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DEVELOPMENT OF DIRECT-INVERSE 3-D METHODS FOR APPLIED TRANSONIC 
AERODYNAMIC WING DESIGN AND ANALYSIS 

I . Introduct ion 

This report summar i zes the activities and accomp 1 i - ; ! ;;;en ts 
associated with Texas A&M Research Foundation Project 5373 which was 
■funded as NASA Grant NAG-1-619 -from the NASA Langley Research Center. 
The project was awarded October 15, 1985 and actively continued until 
August 31, 1989. The primary objective o-f this e-f-fort was the 
development o-f a three dimensional direct-inverse transonic wing design 
and analysis code based upon the TAWFIVE analysis code (Re-f. 1). 

Because o-f its complex nature and the desire to establish proo-f o-f 
concept prior to -final code development, the project was divided into 
two phases. The -first phase developed an inviscid design code, 
established the validity o-f the method, and demonstrated the versatility 
o-f the approach by designing entire wings and discontinuous sections o-f 
wings. The second phase extended the method to include viscous 
interaction and investigated the limits and utility o-f the method. In 
addition, it indicated that it is -feasible to success-ful 1 y design a 
region of a wing which begins aft of the leading edge and which 
terminates prior to the trailing edge. 

1 1 . Personnel 

While the project was officially awarded in October 1985, the 
fiscal paperwork was not completed for several months and the actual 
work was not permitted to start until January 1986. At that time Mr. 
Thomas A. Gaily was assigned to the project as a graduate research 
assistant (GRA). Mr. Gaily remained with the project on a GRA basis 
thru August 1987, and since then he has assisted the project whenever 
needed. In June 1987 Mr. Robert R. Ratcliff joined the project as a GRA 
to conduct the second phase of research. Mr. Ratcliff remained with the 
project thru August 1989. Both Mr. Gaily and Mr. Ratcliff used their 
research work on the project as the basis for their master's theses. 
Mr. Gaily recevied his M.Sc. degree in May 1987, while Mr. Ratcliff 
received his M. Sc. degree in August 1989. 

The principal investigator for this project has been Dr. Leland A. 
Carlson, Professor of Aerospace Engineering. Originally, this entire 
project was to last two to three years, with each phase requiring about 
half of the total time. However, due to the discovery of a spanwise 
oscillation problem during the second phase of the project, the latter 
portion has taken considerably longer than anticipated. The principal 
investigator apologizes to NASA for this delay. 
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in the following seven publications: 
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The research established that: 


<1) The method could obtain invsicid wing designs in v both 
subscritical and supercritical -flow. 

(2) The method could be used to design entire wings or 
noncontiguous regions of the wing on both the upper and lower 
surf aces . 

In addition, it was shown that the method could handle twist, could be 
used to change a wing -from supercr i t i cal to subcritical, and could be 
used to make large surface changes to the original wing. 

In the second phase of the project, viscous interaction was added 
to the design method. In addition, extensive studies of the method and 
comparisons with other codes were conducted in order to verify the 
method. These other codes contained a mixture of similar and different 
coordinate systems, flow solvers, and design methods. Based upon these 
studies, <Ref. 6 ) , it was concluded that the present method and code was 
reliable and accurate. It addition, it was determined that inverse 
methods using similar coordinate systems and flow solvers will yield the 
same wing designs, and that inverse methods having different coordinate 
systems and fuselage represen tat i ons but similar design procedures will 
yield different section profiles. However, the pressure di str i btu i ons 
and lift coefficients in the latter case will be in reasonable 

agreement . 

In addition, extensive studies were conducted to determine the 
approximate limits on wing aspect ratio and leading edge sweep angle 
required for a successful design. Also, studies showing the effects on 
the final design of spanwise grid skewness, grid refinement, viscous 
interaction, the initial airfoil section , and Mach nubmer pressure 
distribution compatibility were conducted. It was determined that: 

(1) Designing at every other spanwise station is the most 
efficient approach . 

(2) A smoothly varying grid is needed at the wing tip for 

accurate design. 

(3) The final designed airfoil sections are independent of the 
initial sections if the direct-inverse junction is moved 
close to the leading edge. 

(4) Boundary layer displacement thicknesses must be included in 

the design process. Otherwise, the designed wing will have 
less lift and different pressure distributions than desired. 

<5) For the conditions considered, wake curvature and 
displacement effects have very little effect on the designed 
airfoil shapes or on the wing pressure distributions. 

(6) Presently, the design of only high and medium aspect ratio 

wings is possible with this code, although preliminary 
approximate results can be obtained for highly swept low 
aspect ratio wings. 
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INVJSC1D TRANSONIC WING DESIGN USING INVERSE 
IN CURVILINEAR COORDINATES 


Thomas A. Gaily* 
Texas A&M University 
College Station, Texas 


Leland A. Carlson** 
Texas A&M University 
College Station, Texas 


ABSTRACT. 

An inverse wing design method has been developed 
around an existing transonic wing analysis code. e 
anginal analyst code. TAWFIVt. has as us core the 
numerical potential flow solver. FLOiO developed by 
Jameson and Caughey. Features of the analysis code 
include a finite-volume formulation: wing and fuselage 

fitted curvilinear grid mesh; and a viscous boundary 
lover correction that also accounts for viscous wake 
thickness and curvature . The development of the inverse 
methods as an extension of previous methods existing for 
design in Cartesian coordinates is presented. Results arc 
shown for inviscid wing design cases m super-critical 
flow regimes. The lest cases selected also demonstrate 
the versatility of the design method in designing an 
entire wing or discontinuous sections of a wing. 

NOMENCLATURE 


h 

H 

J 

Moo 

<loo 

Q 

U,v,w 


- Coefficient of pressure 

- Jacobian of coordinate transformation 

- Jacobian matrix 

- Transpose of inverse Jacobian matrix 

- Freestream Mach number 

- Magnitude of freestream velocity 
Magnitude of local velocity 

- Components of physical velocity vector 

<• .. : . \>£ 


U V W * UOlupuIlCUio vi K'* - 

u\v,W- Components of contravariant velocity vector 

- Angle of attack 

- Ratio of specific heats 

- Differential operator 

- Displacement thickness 

- Displacement thickness due to relofting 

- Trailing edge thickness 

- User specified trailing edge thickness 

- Averaging operator 

- Density 

- Reduced/perturbation potential function 

($> - 4 > + x cos(cr) + y sin(a) ) 

- Potential function 
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INTRODUCTION 


In recent years the importance of transonic flight to 
both military and commercial aircraft and the develop- 
ment of specialized transonic wings for several flight 
research experiments have prompted significant efforts 
develop accurate and reliable computational methods lor 
the analysis and design of transonic wings. Many methods 
of solution have been developed, but among those which 
have shown promise due to their computational efficiency 
and engineering accuracy have been those based upon the 
full potential flow equations in either their conserva i 
or non-conservative form‘d. The TAWF1VE 4 FORTRAN 


• Graduate research assistant. 

•• Professor of Aerospace Engineering, 
Fellow of AlAA ' 


code in particular has proven .o be an exccllen, and 
reliable analysis tool. This analysis code is based upo ^ 
FLO30 finite volume potential flow metho 
developed bv Jameson and Caughey 3 . Among the fea- 
tures of FLO30 are its fully conservative formulation and 
its three-dimensional curvilinear grid. The latter can 
fit around any general combination of fuselage shape 
wing pbnform. 

The purpose of the research described in this paper 
has been to develop a wing design method that is based 
on the existing TAWF1VE analysis code and is compatible 
with the existing computational methods an ^ P r °® r3 Tj 
structure of that code. Of the many wing and a ‘ rr °' 
design methods available 5 -*, the inverse method a* 
developed by Carlson 9-12 was selected for use The 
current work extends the previously developed design 
methods developed for orthogonal gods u . the more 
generalized curvilinear grid system of TA * 

also providing greater design flexibility and versaI ' 1 • 
engineering applications. These last goals were achieved 
by the inclusion of user options for designing either the 
entire wing or only discontinuous wing segments 
shown III Figure 1. The availability of this option is 
useful to engineers who are typically faced '*> lh des, 8 - 
ning around regions where the wing geometry ma>be 
fixed by constraints other than aerodynamic consider- 

ations. 




Associate 




wing ANALYSIS Methods 

Potential Flqw «fr|. Tf 

The inviscid poiemial analysis of TAWF1VF i< 

ZlT" ‘Vtf »’»«“" FLO30 df .elofwd by CiufK»y 
FLO10 •„ . ,f« < “"»'=« Oe«rip,io„ S*S 
, “ code and '■* theoretical basis the reader it 

referred to Caughey and Jameson's papers and some 
earlier developmental work by Jameson ,4 ‘ 15 a hri^r 

and Cr, to ,,0 D n "• , presenIed here 10 Provide for completeness 
nnd to prov.de a background for the inverse desier 
developments which will be discussed in detail. 

live form 30 w S h O VK eS T P ° temia ' equa,ion in conserve 
e form which when transformed from Cartesian conr- 

mates to generalized curvilinear coordinates is. 

0>hV) ( ♦ (phV), ♦ <phW) f - 0 (j) 

where the subscripts denote differentiation with respect to 

velochfes 'are ' related*"* *** £' a " d f ' The c °rr tra v ar .an t 
d , • tl\ r ? d the phvsical velocities and the 

derivatives of the potential function by: 


V> 


H’ 1 {•"]- (»%,-! {Jj 

and H is the transformation matrix defined by: 
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x< 


y £ y n 
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( 2 ) 


( 3 ) 
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. ThC n “ m . ericaI approach used in FLO 30 is a finite 

. Ume t . echniQue \ To understand this approach, consider 

svstem^n lW ° d * menSl0naI case represented by the grid 
system shown in Figure 2. J h 
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i+l 


Finite- Volume Cell Location 

Smen a t Sh unde Ube ' he figUre indica,es ** area 

sid^Tb «n h C ° nS,dera,,on - The flux of fluid through 
n * * • be approxima ted by the average of the 

c d The^'n and b .. Wi,h Simibr r « ulK fdr ** side 

area centered a . n d‘ n ' be X . direction for the elemental 
area centered at grid point i,j js then: 

(phU)^ » f(phU a + phU b ) - (phU c + phU d )] / 2A( 

or in the notation of Caughey and Jameson, 

* 

W>U) f - »n&((pV) 


* Veragin « ,nd * indicates differentiation 
,"„* b ' n ' nd,q *'* d «*'rections which are defined as follows 
(allowing ■ | ); 

- Ui-U.k) 

tu m. “ <U 'H.j.k * Uj.j j k )/; 

(l ij.k ” + ^i-J.j+j. 


. etc. 


+ U i-i,j-i.kV4 


averir^nV 1 ex,ended 10 •be other flu* components and to 

numeral 0V " SUrfaces in Ihree dimensions, the 

numerical potential equation is of the form: 

V n ^ ( (phU) ♦ Pj-rtyphV) * p f(; « f (phW) - o 

the r T !/’ nd , ,be n “ X quamilies Phv, and phW at 

for fhe , H Ume ver,ices (i e - po,ms a - b. c, and d 
EouatV °„ d,m ! nS, ° nal ca5e >- ’« “ necessarv t0 evaluate 
Equations (2) through (4). The derivatives in these 
expressions can be expanded by the same volume avcrac- 
mg approach used above, thus: 5 




W) 

"ffty*) 

P<l7*r<*> 


x ( = p-3f¥ x ) 
y ( = 

z f “ "ij&w 

The h ‘ ermS - f ° r ,be °‘ her information metrics. 

The above expressions, being centered at grid midpoints 
will involve the values of the potential and grid position 
at grid points which are known from the previous poten- 
tial solution and the grid geometry, respectively. 

When solving transonic flows it is necessarv to 

upstream"^* Sqlu,,on . al 8° ri,hm «>me form of superionic 
upstream dependence in order to account for both the 

Physical nature of the flow and the viscous nature of 

durett WaVeS ’ reSpeCI,vely ' Caughey and Jameson intro- 
duced upwnuimg by , he addition of terms into their 

When the "no en ”‘ equa,i0n which are °"Iy non-zero 
hen the flow is supersonic. Also, the finite volume 

fie C .d m slf Xhib l tS a ,endenCV f ° r uncoupling of the now 
r f d , S °'“ npn Pe,ween alternating grid points. As a 
„ . addmonal «rms are included in the numerical 

potential equation. The final numerical equation which is 

has Ve the for^ L03 ° WbCn ' h “ e " rms have been includ c d 
PqfS&hU+P) + Pfff|j(phV+Q) + p, .f f (phW+R) 

■ + + * ^Q^/2) - 0 

where P, Q. and R are the upwinding terms and Q, n , 

30 are the decoupling terms. 

Compu^tional Grid Geom^ffy 

n ,J he computational grid used by FLO30 is a body 
1TT ,#M !'. curvilinear mesh constructed about 
a w,ng fuselage combination. The number of grid points 
omposing the computational domain is typically 40 x 6 x 
n Xl6 ' -° r 160 * 24 x 32 for the number off 

’’ a " d 1 po,n ll ,n ,he coarse, medium, and fine grids, 
respectively. The grid is conformally mapped to the 
wing and fuselage surfaces as can be seen from the p,o! 
of surface grid lines shown in Figure 3. 

The grid is formed around spanwise airfoil sections 
in a similar manner in which *C grids are mapped to 
airfoils in two-dimensional analysis. In addition each 
spanw.se computational plane is also conformally wrapped 
about the fuselage surface and a line extending forward 
from the fuselage nose. 




Figure 3. Surface Grid Point Geometry 

A final set of grid surfaces are generated beneath 
the wing and fuselage surfaces and beyond the symmetric 
plane in order to aid in the formulation of both the 
finite-volume numerical flow equations and the flow 
tangency boundary conditions upon these boundaries. 
The grid points composing the "ghost*' surfaces are 
calculated from linear extrapolations of the computation 
grid lines from inside the physical domain. 

Boundary Conditions 

Since the governing potential equations are written in 
terms of perturbations from free-stream conditions, the 
subsonic, far-field requirement that the flow return to 
the free-stream velocity and direction is satisfied by 
setting the perturbation potential equal to zero on the 
side and upstream boundaries. The downstream boundary 
condition is a "zero" order extrapolation of the potential 
(constant potential assumption) to the outflow boundaries. 

A flow tangency condition is applied along both the 
wing and fuselage solid surfaces by setting the normal 
contravariant component of the velocity vector to zero on 
the surfaces. This condition provides an equation which 
when approximated by a finite-difference expansion 
about the surface grid points can be used to set a value 
for the perturbation potential on the "ghost" grid points 
below each surface. Note that this finite-difference 
boundary condition differs in formulation from the 
finite-volume solution algorithm of the governing 
equations. As a result, it is possible to impose flow * 
tangency using the finite-difference technique yet still 
have a slight normal surface velocity when performing 
the finite-volume calculations. Since it is essential to 
have accurate boundary conditions at the wing surface in 
order to generate accurate solutions, a second condition is 
imposed upon the wing surface. This additional condi- 
tion involves reflecting the flux quantities calculated by 
the flow solver for the cell centers directly above the 
wing surface to the "ghost" cell centers beneath. The 
reflected normal fluxes then cancel each other out in the 
residual expression and a net zero flow is obtained 
through the surface. Similarly, a zero flux condition is 
applied at the half -body symmetric plane, limiting 
solutions to symmetric, non-sideslip cases. 

The trailing edge slit boundary is not an actual limit 
to the physical domain as the other boundaries are, but is 
simply an artificial boundary created by unwrapping the 
physical plane into the computational domain. The only 


conditions which need to be imposed at the slit is that 
the flow velocities, and thus pressure, be continuous 
across the cut. The flow potential, however, will have a 
discontinuous jump across the wake which is proportional 
to the sectional wing lift coefficient. 

INVERSE WING DESIGN METHODS 

As stated previously, a direct-inverse approach to 
wing design was selected for incorporation into the 
TAWFIVE code. The direct-inverse method derives it* 
name from the division of the design wing surface into a 
fixed geometry leading edge region, where flow tangency 
boundary conditions are imposed, and an aft, variable 
geometry section where pressure boundary conditions are 
enforced. The pressure boundary where the user speci- 
fied pressure distributions are imposed does not extend 
forward to the leading edge due to difficulties of 
enforcing this type boundary condition near the beginning 
of an airfoil section. This restriction on the size of the 
pressure specification region does not seriously reduce the 
versatility of the design method since the leading edge 
regions for most airfoils are similar, and it is relatively 
easy to select a leading edge geometry which will 
produce the desired Mach number or pressure values at 
the beginning of the inverse region. In addition, specific 
leading edge shapes may be required due to other design 
constraints such as the necessity to house a leading edge 
flap or slot system. 

Pressure Boundary Condition 

In the inverse design regions on the wing, a pressure 
boundary condition will be specified rather than the flow 
tangency condition used in analysis zones. In formulating 
this boundary condition it is necessary to relate the user 
specified pressure coefficient, C p , to the current 
perturbation potentials at inverse design grid points. 
Consider the full potential equation for the pressure 
coefficient: 


. i 2 i[i Wood 

1Moc [L 



where: Q 2 


If it is assumed that the pressure coefficient is 
primarily a function of the chordwise component of the 
velocity, u, and only slightly affected by the vertical and 
spanwise components of velocity, v and w, then a stable 
approximation is made by time lagging the latter two 
velocities in the boundary condition expression. This 
assumption is true everywhere except near the leading 
edge; but since the inverse design boundaries have 
already been restricted to regions well behind the leading 
edge, the simplification is justified. The value of the 
local velocity, u, can then be calculated from the above 
expression in terms of the desired pressure coefficient 
and the current values for the vertical and spanwise 
velocities. In addition, the velocity u can also be 
calculated from the perturbation potentials using the 
relations of Eq. (2). Defining Jy to be the elements of 
the inverse transpose of the Jacobian matrix, H, the two 
equations for u yield: 


*ll^ +J 12V- J 13*f 


t 1 

7 


_L_2|7 7MooC p \ 

‘h-DMoo^l +— 2 )- I 

< * ilf ♦ (S ') 2 


cos(a) (5) 


already bJn , P j * nd verncal now velocities have 
condV , “ un ' ed W h* cons,an < in the boundary 
condition relation, it is consistent to make the same 

approximation in the above expression with respect to the 
span wise and vertical derivative terms. i„ andT This 
assumpnon ,s similar to the previous one.’and leaSs to an 
explicit expression for the potential at one point. 

™. e rin !. ,e difference approximation used involves 
“CS.7 i! * *bout° the 

central rti rLl.' • , ( der,va,,ve >s determined bv a 

centra) difference involving , he preceding and following 

fh mT V 1 ' The ’ a " d f d '"vatives are found at 
nrerert d ' P0,T '} by av « r aging the derivatives from the 
nn^ r, and f0,,0win8 grid P° inIS round by a thee 
respec,ive C | Wa FL a f / e K ,f:i1 differen « ^ima.ions! 
Pressure IpecifSn S % tS 5 ? The 

^maZT^' eXPreSSi ° n ° b,ained Wilh * b '” ^nite 


. /a"* 1 " 

J ll^i.j.k - *i-lj fk ) 



+ J i2[w5i + - Wij“ , k * 

, n n 

* *ij-2,k + ^i-l,j-2,kj /< 

J J3^i,j.k+] + <*i-l,j,kvl * ^i.j.k-I - ^i-l.j.k-J 

• F (CPi-4, k ) 

value^of ^ SUpereCri P“ n and refer to current 

values of the potential and the new values of the 

reTpZveh i. b >’ ' he boundary condn.on 
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• r a Q ^ evaluated using the pressure coefficient 
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The potential values at n+1 in the direct region are 

inverae h” 1 * h"* S ' nCe ‘ hey d ° n0t change when 'he 
inverse boundary condition is applied; i.e. <* n+1 . ^ All 

calculated" and* °" '*1 ' nVerSe boundar y can then be 
“e smln I*,, •* ,h * Spanwise and vertical derivatives 

X!; J, Primarily be functions of the pressure 
coefficient at grid point i-i and the value of the 
potential at grid point i-1. 


Mfm! h !,-I? nly c . oncern with usin 8 'his mid-point specifi- 
cation scheme is that the current method of calculating 

entered S dlffe a * a "“T fr0m FL ° 30 “*« a «" d 3 
Thi! dfffe enCe * Cheme for the stream wise derivative. 

soeciffed Z Ce , C ° U J d POtemia,ly aIlow a Pressure to be 
P f5 ,ed CO " ectly *>ut still have a significantly different 

lion Methods f u° m FL ° 3 ° dU * 10 ,he inconsistent calcula- 
tion methods. However, as shown on Figure 5, where the 

compared f ° f 3 typi «' ^w solution are 

thTnA, d Ki f0r C tW ° different calculation techniques, 
this possible error has nert been significant in practice. 


Figure 4. Point Dependancc and Location 



Figure 5. Comparison of Pressure Calculation Methods 
Surface Calculations 


me inverse boundary conditions drive the flow 
field to a converged solution, it is necessarv l0 
periodically calculate the location of the new displacement 
surface and to regenerate the computational grid about 
this new geometry so that the pressure boundary surface 
will correspond to the physical boundary surface. Each 
new surface can be found relative to the* previous surface 
rom an integration of the wing surface slopes. However 
the surface slopes must first be calculated from the 
current flow field solution using the flow tangency 
boundary condition which in curvilinear coordinates is: 


VT 


x VF 


where V is the contravariant velocity vector and VF is 
the gradient of the surface function with respect to the 
curvilinear coordinates. Note this condition is a direct 
analog to the same condition expressed in physical space. 

A more useful expression can be obtained bv 
expanding the above equation to: 



expression can be solved for the new chordwise 


X \X 



(SO , 


cotwtant* nT^he * curren^slopes on the ^ r '* su ^' equal 
zero and a simplified How tangency condition results. 

V 

u 

r wing 

The above expression has been applied th J e) e ™; 
pu.ztional ^- r ace plane .n order ^to ^ #pproach 

true when appl»e stable iterative surface 

astss.’ss — — 

surface. 

Tn calculate the relative surface slopes, it « 
necessary to accurately determine the 

C nrbTThe V tXofVeS/e« al.'^ a simple finite 
difference calculation of * e “ «' , 
iT^accurate^me'thod ^“Implemented which uses the 

under the assumption that the resid FLO30 can 

surface points. The residual expression from FLO30 
be written in finite volume form as. 

n^OthU) ♦ #f(ty*V) ♦ ,^W) term$) _ 0 
The "other terms' in the above expression involve the 

rtfsur- ... 

following development. 

The desired velocities can also be written in 
finite volume form as: 

V - phV - P W (/>hV) and U - phU - #»<i,fO> hU > 

By simple manipulations, the normal velocity can be 
obtained from the residual expression as: 

2„ w w.v) - <« 

where the subscript ,-l refers to the values at grid cell 
centers above the wing surface. 

In order to use Eq. (6) to find the desired surface 
velocity ratio, it is necessary to know^ ^ 

velocity components at the ghost . manner 

wi „ g surface. Thw values -n be obta.ned m ^ ^ 

consistent w.th by speedy, g points 

values to equal the values at cur . of 

j^SScy^ 

th?ch the calculated slopes should of course be zero. 

With the contravariant velocities known, an integra- 
tion of ^(6) through the inverse 

the leading edge to the trailing * d8 * J wing surface 

?TS5J^ Ai V-VSiT 

EXDrcssed as changes in the computational coordinate fl, 
”d a e converted to surface displacements tn the 
physical plane via the, local grid tmnJorm«.on The 

physical plane displacements are coincident 
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produce unpredictable problem* « the gr.o 
flow calculation portions of FLU3U. 
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boundary condition for the 1 viewed as the proper 

assuring trailing edge closure procedure for 

selection of a lending edge *»** * region in order 
systematically mod ‘^"V«i£g edge thickness is called 
to achieve some des red 8 8 h been incorpor- 

relofting. Such a relof, ting h “ order t0 both 

ated into the Present edge crossover and to 

prevent the problems of tra.li g * trailing edge 

allow the user the opt.on of specifying ^ design 
thickness as an additional d«jg . , applications 

feature should a leading edge 

SZ*£STZ ToZ^l have to be performed by 
the user. 

Two methods of 

The first method is * * ,mp ( * e hclp of Figure 7. The 
method can be v,s “ al, “ origm al leading edge geometry 
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calculated for the tnv ace haJ a tra iling edge 
modification, this n ^ A wcre specified by 

thickness of J f * f S^kf hfe to be relof.ed or 
changed^ ^"the present scheme, in order to obtain 





Figure 7. Relofting to Force Trailing Edge Closure 

desired thickness, a displacement thickness, 6 r , is added 
to the current design surface. This thickness has a 
distribution from the leading to the trailing edge and is 
determined by the formula: 

M*) - < A t “ A) < x/c) 

where c is the chord length of the local airfoil section. 
The total displacement for a surface update is then the 
sum of the two displacements, 6(x) and f r (x). When 
both the upper and lower surfaces are designed simulta- 
neously, the displacement magnitudes determined by 
relofting are divided between the two surfaces so that 
half is added to the lower surface and half to the upper 
surface. 

The second relofting method uses the same approach 
as the first for the aft inverse regions, but modifies the 
leading edge region by a proportional thining or thicken- 
ing of the surface ordinates. This approach can be 
expressed by: 

y n+! (x) « yj n+1 y n (x) / yj n 

where the j subscript refers to the ordinate at the 
direct-inverse junction determined from the linear 
relofting of the aft regions. Note that this method will 
produce leading edges in the same family of shapes and, 
for example, allow the design from a NACA 0012 airfoil 
to a NACA 0006 airfoil (see Test Case F). 

RESULTS 

A variety of different test cases were run as 
verification of the current design method. These cases 
involved both subcritical design and supercritical design 
over section geometries selected to test the versatility of 
the input and design control logic. In this section results 
from three of the more significant test cases wilt be 
presented. The results shown were obtained on a 
medium grid having 81 streamwise, 13 vertical, and 19 
spanwise points with 1 1 spanwise stations and 53 points 
on the wing at each station; and in all cases the 
maximum change in the reduced potential was reduced at 
least three orders of magnitude. Thus, the results do not 
represent ultimate convergence but should be represent- 
ative of "engineering accuracy". 

The planform selected for the test cases was the 
Lockheed Wing A wing-body. The wing for this config- 
uration has a quarter chord sweep of 25 deg., a linear 
twist distribution ranging from 2.28 deg. at the wing 
body junction to -2.04 deg. at the wing tip, an aspect 
ratio of eight, and a taper ratio of 0.4. The last two 
values are based upon the wing without fuselage. 
However, instead of the supercritical sections normally 
associated with Wing A, the initial airfoil sections at each 
span station were assumed to be composed of symmetric 
NACA four digit airfoil sections. 

The target pressure distributions used in the design 


regions for the first two test cases were selected to yield 
airfoil shapes thicker in the aft portions of each section; 
and, at supercritical conditions, to yield on the upper 
surface weaker and more forward shock waves than those 
which would normally occur on a NACA 0012 section. 
On the lower surface, the target pressure distributions 
were selected to have either a favorable pressure gradient 
or fairly constant pressure plateau over much of the 
lower surface. 

For the last test case, the pressure distribution was 
obtained from analysis solutions of an assumed wing 
geometry. The intent of this cases is to verify the 
relofting procedures and show the ability of the current 
method to make large surface changes in going from a 
thick wing to a thin wing (approximately 12 percent to 6 
percent thick respectively). 

All cases were for a freestream Mach number of 0.8 
and an angle of attack of two degrees. In each case, the 
pressure distribution was specified in the design regions 
from the 15% local chord location to the trailing edge 
and used as the boundary condition in these inverse 
regions starting with the first iteration. Normally, three 
hundred SLOR iterations were executed prior to the first 
design surface update calculation; and subsequently, 
surface updates were computed every fifty cycles. 
Usually, the solution was considered converged and 
terminated after 550 total iterations for the first two 
cases and, due to the large amount of relofting required, 
after 950 iterations for the last case. 

Test Case C 

The inverse design regions for Case C. which was an 
attempt to design both upper and lower surfaces on two 
noncontiguous regions of the wing at supercritical 
conditions, are shown on Figure 8; and a comparison 
between the initial pressure distribution associated with 
NACA 0012 sections and the target pressures for two 
sections is portrayed on Figure 9. As can be seen, the 
target pressure distribution essentially eliminates at 
inboard stations the upper surface shock wave present on 
the original wing; and at outboard stations it weakens the 
shock and moves it forward. In addition, significant 
changes in the lower surface pressure gradients are 
evident. Also shown on Figure 9 are the pressures 
computed by the program at the end of the inverse 
design procedure (denoted as "design pressures"). These 
pressures are in excellent agreement with the target 
pressures, which indicates that the method is satisfying 
properly the desired inverse boundary conditions. 

The corresponding designed airfoil sections for this 
case are shown on Figure 10. Even on the expanded 

scale, the agreement between the designed and target 
surfaces is excellent at all design stations. However 
trailing edge closure was not enforced for this case; and 
there is at the boundary stations some departure between 
the designed surfaces and the target surfaces near the 
trailing edge. It is believed that this slight difference is 
a ramification of the change in spanwise slopes near the 
trailing edge between the direct and inverse regions. 

In any event, the pressure distributions resulting 
from an analysis of the designed surfaces shown in 
Figure 10 are in excellent agreement with the target 
pressures, as can be seen on Figure 11. In addition, the 
section lift coefficients at the various design stations are 
in very good agreement with the target coefficients. 
Based upon these results it is believed that the present 
method can adequately design/modify nonadjacent regions 
of a wing in transonic flow. 
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However, as depicted on Figure 14, only the design 
surfaces change form the original shape; and these 
surfaces are in reasonable agreement with the target 
profiles. 

Finally, Figure 15 compares analysis results obtained 
for the designed wing with the target pressures. Even 
for this complicated case, the agreement between the two 
distributions and between the actual and target lift 
coefficients is excellent. 


inverse design regions 
DESIGN CASE E 



Figure 12. Design Case E 



Figure 13. Comparison of Initial Pressures with 
Targei Values (Case E) 



Figure 13. Continued 



Figure 14. Comparison of Designed Sections with 
Original and Target Sections (Case E) 
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TEST case f 

The final test case was selected to demonstrate the 

ss « f r °*s sss ™ rs. 

large changes in the leading edge reg.on through 

*s a.’sj'Ss 

case are shown in Figure 16 where the wing «h'ckne« 
varied from 12% to 6% between the wing too and 80% 
span location and was constant going outward to the tip. 
The input design pressures were for a constant 6% thick 
wing. 


The first attempts at this design used the linear 

leading edge relofting Procedure e "?J design surfaces 
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Figure 17. Comparison of Initial Pressures with 
8 Target Values (Case F) 
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Figure 17. Continued 


Figure 18, Continued 
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Figure 18. Comparison of Designed Sections with 
Original and Target Sections (Case F) 


Figure 19. Comparison of Pressures from Analysis of 
Designed Wing Target Distributions (Case F) 
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CONCLUSIONS AND SUGG ESTIONS 
FQR FUTURE WORK 

A direct-inverse wing design method has been suc- 
cessfully incorporated into the TAWFIVE transonic 
wing-body analysis computer code. The resultant code is 
capable of designing or modifying wings at both tran- 
sonic and subsonic conditions and includes the effects of 
wing-body interactions. A series of test cases have been 
presented which demonstrate the accuracy and versatility 
of this inverse method. 


Inclusion of viscous effects via the addition of the 
wing surface displacement thickness and wake thickness 
when performing wing design has been accomplished but 
not completely verified. Additional work will be 
required to run a sufficient sampling of test cases for 
evaluation of this design mode. The unique problems 
associated with viscous design and the ef / cc ” ° r * * 
various viscous correction models available in TAWFI t 
would be the subject of a continuing research effort. 

The development and evaluation of alternate methods 
of surface relofting are also topics for which continued 
research is suggested. The current method of relofting 
restricts the user to a family of leading edge geometries 
which can be constructed by the linear rotation of the 
initial shape. The option of using other relofting 
methods would extend the family of available shapes and 
add versatility to the design method. 
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in Curvilinear Coordinates 
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An Invent wing design -hod bes been developed 

original analysis code, TAWFIVE, bas as its core ■ * j tonntt i»tlon, an SLOR soluUon 

Caughey and Jameson. Features of tbe an y development of the Inverse method as an ex- 

scheme, and a wing and fuselage filled, curv near g ^oordinetes b ©reseated Result* are shown for in- 

tension of previous methods existing for design I. tb e versaUlity of 

Viscid wing design cases In supercriticd flow regimes. Tbe test 
Ibc design method in designing an entire wing or discontinuous sections of a wing. 
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Nomenclature 

= chord length 
= coefficient of pressure 
r Jacobian of coordinate transformation 
= Jacobian matrix 

= transpose of inverse Jacobian matrix 
= freest ream Mach number 
= upwinding terms 
= decoupling terms 
= magnitude of local velocity vector 
= magnitude of freestream velocity vector 
= wing surface function 
= components of physical velocity vector 
= components of contravariant velocity vector 
= contravariant velocity vector 
= Cartesian coordinate directions 
= angle of attack 
s ratio of specific heats 
= differential operator 
= displacement thickness 
= displacement thickness due to relofting 
= trailing-edge thickness 
= user-specified trailing-edge thickness 
ss decoupling factor 
= averaging operator 
= curvilinear coordinate directions 
= density 

= reduced/perturbation potential function 
= potential function ($ = <#>+* cosa+y sina) 


Introduction 

I N recent years, the importance of transonic flight to both 
military and commercial aircraft and the development of 
specialized transonic wings for several flight research ex- 
periments have prompted significant efforts to develop ac- 
curate and reliable computational methods for the analysis 
and design of transonic wings: Many methods of solution have 
been developed, but among those that have shown promise 
due to their computational efficiency and engineering ac- 
curacy have been those based upon the full, potential flow 
equations in either their conservative or nonconservative 
form. 1 ' 1 The TAWFIVE 4 code in particular has proven to be 
an excellent and reliable analysis tool. This code is based upon 

Presented as Paper 87-2551 at the AIAA 5th Applied Aerodynanucs 
Conference, Monterey, CA, Aug. 17-19, 1987; received Sept. 11. 1987, 
revision received Oct. 27. 1987. Copyright © 1987 American Institute 
of Aeronautics and Astronautics, Inc. All rights reserved. 
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t Professor of Aerospace Engineering. Associate Fellow AIAA. 


the FLO30 finite-volume potential flow method that was 
developed by Caughey and Jameson. 5 Among the f « tur “ ° f 
FLO30 are its fully conservative formulation and its three- 
dimensional curvilinear grid. The latter can be fit * rou "l5’ y 
general combination of fuselage shape and wing planfonm 
The purpose of the research described in this paper has been 
to develop a wing design method that is based on the existing 
TAWFIVE analysis code and is compatible with the existing 
computational methods and program structure of that 
Of the many wing and airfoil design methods available, th 
inverse method as developed by Carlson,* 10 Anderson and 
Carlson 11 and Weed ct al. 12 was selected for use. The current 
work extends the previously developed design methods 
developed for orthogonal grids to the more generalized cur- 
vilinear grid system of TAWFIVE, while also providing 
greater design flexibility and versatility for engineering ap- 
plications. These last goals were achieved by the inclusion of 
user options for designing either the entire wing or only 
discontinuous wing segments as shown m Fig. • 
availability of this option is useful to engineers who are 
typically faced with designing around regions where the wing 
geometry may be fixed by constraints other than aerodynamic 
considerations. 

Wing Analysis Methods 

Potential Flow Solver 

The inviscid potential analysis of TAWFIVE is per- 
formed by the program FLO30 developed by Caughey^nd 
Jameson 3 n For a complete description of the FLO30 code 
and its theoretical basis, the reader is referred to Caughey and 
Jameson’s papers and some earlier dcvel °P ment ^ 

Jameson. 1415 A brief description is presented here to provide 
for completeness and a background for the inverse design 
develonments that will be discussed in detail. 

FLO30 solves the full potential equation in conservative 
form that when transformed from Cartesian coordinates 
generalized curvilinear coordinates is 

( P hU) ( + (phV)i+(phW) f =° 0 ) 

where the subscripts denote differentiation with respect to the 
curvilinear coordinates and f. The contravariant velocity 
are related to the physical velocities and the derivatives of th 
potential function by 


( 2 ) 




r«i 
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= [H T H]-‘ 
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where H is the transformation matrix defined by 

r *t v 

H« y i y, y r (3) 

^ z t - 

with n- IHI. An equation for the local density can be ob- 
tained from isentropic relations as 




1 + ** ■ ■ ] A/1(1 - u 2 — v 2 - Vi* 2 ) 


i'b-0 


(4) 


The numerical approach used in FLO30 is a finite-volume 
technique. To understand this approach* consider the simple 
two-dimensional case represented by the grid system shown in 
Fig. 2. The dashed cube shown in the figure indicates the area 
element under consideration. The flux of fluid through side 
a-b can be approximated by the average of the fluxes at point a 
and b with similar results for the side c-d. The net flux in the x 
direction for the elemental area centered at grid point ij is 
then 



= (f4* vij* vi.* + £4* n j- vi,k 


(phU) i = 


HphU.+phU b )-(phU'+phUt)] 

2A£ 


or, in the notation of Caughcy and Jameson, 

(phW^n^iphU) 


(5) 


(6) 


+ Ut- hj + vi.* + Ut - vi j- vi,*)/4 (7c) 

When extended to the other flux components and to averag- 
ing over cube surfaces in three dimensions, the numerical 
potential equation is of the form 

(phU) +Mf,5 r (p/iH0 =0 (8) 


where 6 and p are differentiating and averaging operators in 
the indicated directions that are defined as follows (allowing 
A£ = Atj = Af = 1): 

— (Ui+ViJ.k ~ Ui-KJ.k) P*) 

(PiO) ij.k — ( ViJ.k + ^4- ViJ.k )/2 (7b) 





To find the flux quantities phi), phV, and phW at the 
finite-volume cell vortices (i.e., points a, b, c, and d for the two- 
dimensionaJ case), it is necessary to evaluate Eqs. (2-4). The 
derivatives in these expressions can be expanded by the same 


volume averaging approach 

used above, thus 




(9a) 



(9b) 



(9c) 


with similar terms for the other transformation metrics. The 
above expressions, being centered at grid midpoints, will in- 
volve the values at grid points of the potential and grid posi- 
tion, which are known from the previous potential solution 
and the grid geometry, respectively. 

When solving transonic flows, it is necessary to include 
some form of supersonic upstream dependence and artificial 
viscosity in the solution algorithm in order to account for the 
physical nature of the flow and the viscous nature of shock 
waves, respectively. Caughey and Jameson introduced 
upwinding by the addition of terms into their potential 
numerical equations that are only nonzero when the flow is 
supersonic. These terms also introduce a numerical error that 
has the form of a viscous term. Additional terms are also in- 
cluded to correct a tendency of the flowfield solution to un- 
couple between alternating grid points. The final numerical 
equation, which is solved by FLO30 when these terms have 
been included, has the form 

(phU+P) ( P hV+ Q) 

+ P ilt 6{(phW+R)-e ( f 

+ r' 2 >=° < 10 > 

where P, Q t and R are the upwinding terms; Q tn , Q^ t Qn> 
and are the terms reducing odd-even decoupling; and i is 
a factor determining the amount of decoupling (typically 
€ = 0.25). 
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Computational Grid Geometry 

The computational grid used by FLO30 is a body-fitted, 
nonorthogonaJ, curvilinear mesh constructed about a 
wing/fuselage combination. The number of grid points com- 
posing the computational domain is typically 40x6x8^ 
80x12x16. or 160x24x32 for the number of £, i>. and f 
points in the coarse, medium, and Fine grids, respectively. The 
grid is conformally mapped to the wing and fuselage surfaces 
as can be seen from the plot of surface grid lines shown 

in Fig. 3. . 

The grid is formed around spanwise airfoil sections in a 
similar manner in which ”C” grids are mapped to airfoils in 
two-dimensional analysis. In addition, each spanwise com- 
putational plane is also conformally wrapped about the 
fuselage surface and a line extending forward from the 
fuselage nose. The reader is referred to Refs. 3 and 13 for ad- 
ditional details on the method of grid generation. 

An additional set of grid surfaces are generated beneath the 
wing and fuselage surfaces and beyond the symmetric plane in 
order to aid in the formulation of both the finite-volume 
numerical flow equations and the flow tangency boundary 
conditions on these boundaries. The grid points composing 
the “ghost” surfaces are calculated from linear extrapolations 
of the computation grid lines from inside the physical domain. 

Boundary Conditions 

Since the governing potential equations are written in terms 
of perturbations from freestream conditions, the subsonic, 
far-field requirement that the flow return to the freestream 
velocity and direction is satisfied by setting the perturbation 
potential equal to zero on the side and upstream boundaries. 
The downstream boundary condition is a “zero “-order ex- 
trapolation of the potential (constant potential assumption) to 
the outflow boundaries. 

A flow tangency condition is applied along both the wing 
and fuselage solid surfaces by setting the normal contravanam 
component of the velocity vector to zero on the surfaces. This 
condition provides an equation that, when approximated by a 
finite-difference expansion about the surface grid points, can 
be used to set a value for the perturbation potential on the 
“ghost” grid points below each surface. Note that this finite- 
difference boundary condition differs in formulation from the 
finite-volume solution algorithm of the governing equations. 
As a result, it would be possible to impose flow tangency using 
the finite-difference technique yet still have a slight normal 
surface velocity when performing the finite-volume calcula- 
tions. Since it is essential to have accurate boundary condi- 
tions at the wing surface in order to generate accurate solu- 
tions, a second condition is imposed on the wing surface. This 
additional condition involves reflecting the flux quantities 
calculated by the flow solver for the cell centers directly above 
the wing surface to the “ghost” cell centers beneath. e 
reflected normal fluxes then cancel each other out in the 
residual expression and a net zero flow is obtained through the 
surface. Similarly, a zero flux condition is applied at the half- 
body symmetric plane, limiting solutions to symmetric, non- 
sideslip cases. _ A 

The trailing-edge slit boundary separating the upper ana 
lower half planes is not an actual limit to the physical domain 
as the other boundaries are, but is simply an artificial bound- 
ary created by unwrapping the physical plane into the com- 
putational domain. The only conditions that need to be im- 
posed at the slit are that the flow velocities, and thus pressure, 
be continuous across the cut. The flow potential, however, 
will have a discontinuous jump across the wake that is propor- 
tional to the sectional wing lift coefficient. 



Inverse Wing Design Methods 
As stated previously, a direct-inverse approach to wing 
design was selected for incorporation into the TAWFIVE 
code. The direct -inverse method derives its name from the 
division of the design wing surface into a fixed geometry 
leading-edge region, where flow tangency boundary condi- 
tions are imposed, and an aft, variable geometry section where 
pressure boundary conditions are enforced. The pressure 
boundary where the user-specified pressure distributions are 
imposed docs not extend forward to the leading edge due to 
difficulties of enforcing this type of boundary condition near 
the beginning of an airfoil section. This restriction on the size 
of the pressure specification region does not seriously reduce 
the versatility of the design method since the direct region can 
be fairly small (as little as 3*7© chord), the leading-edge regions 
for most airfoils are geometrically similar, and it is relatively 
easy to select a leading-edge geometry that will produce the 
desired Mach number or pressure values at the beginning of 
the inverse region. In addition, specific leading-edge shapes 
may be required due to other design constraints such as the 
necessity to house a leading-edge flap or slat system. 


Pressure Boundary Condition 

In the inverse design regions on the wing, a pressure bound- 
ary condition will be specified rather than the flow tangency 
condition used in analysis zones. In formulating this boundary 
condition, it is necessary to relate the user-specified pressure 
coefficient C p , to the current perturbation potentials at in- 
verse design grid points. Consider the full potential equation 
for the pressure coefficient: 



f[-^- 




(id 


where q 1 = u 7 + u 2 + w 2 . 

If it is assumed that the pressure coefficient is primarily a 
function of the chordwise component of the velocity u and 
only slightly affected by the vertical and spanwise components 
of velocity v and w, then a stable approximation is made by 
time lagging the latter two velocities in the boundary condition 
expression. This assumption is true everywhere except near the 
leading edge; but since the inverse design boundaries have 
already been restricted to regions behind the leading edge, the 
simplification is justified. The value of the local velocity u can 
then be calculated from the above expession in terms of the 
desired pressure coefficient and the current values for the ver- 
tical and spanwise velocities. In addition, the velocity u can 
also be calculated from the perturbtion potentials using the 
relations of Eq. (2). Defining J v to be the elements of the in- 
verse transpose of the Jacobian matrix H, the two equations 
for u yield: 



*3 


Jn 4>i + - 


— cosa 


( 12 ) 
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Since the span wise and vertical flow velocities have already 
been assumed to be constant in the boundary condition rela- 
tion, it is consistent to make the same approximation in the 
above expression with respect to the span wise and vertical 
derivative terms 4, and This assumption is similar to the 
previous one, and leads to an explicit expression for the poten- 
tial at one point. 

The finite-difference approximation used involves expand- 
ing the derivatives of the potential about the midpoint /— Yi t j, 
k. The i derivative is determined by a central difference in- 
volving the preceding and following grid point values. The t) 
and f derivatives are found at the midpoint by averaging the 
derivatives from the preceding and following grid points 
found by three-point-backwards and central-difference ap- 
proximations, respectively. Figure 4 shows the point 
dependence and pressure specification point for this method. 
The resulting numerical expression obtained with these finite 
approximations is: 

+ J,1 [3 (*&' + u.k ) M + V-u- U ) 




= F(C p .,_ WJt ) (13) 

Here,the superscripts n and n + 1 refer to current values of the 
potential and the new values of the potential being imposed by 
the boundary' condition, respectively. Also, the term 
F{C Pii _ njc) is the right-hand side of Eq. (12) evaluated using 
the pressure coefficient specified at point / — Vi, k. Solving 
the above expression for the potential at point ij t k yields 


4> 


U.k 


1 

J n + V*J i2 




/— 1 J.k 






“ ^13 ( + 1 + *?- IJ.k + 1 4>i, J.k - 1 


The potential values at n + 1 in the direct region are known 
initially since they do not change when the inverse boundary 
condition is applied; i.e., + 1 =<f>”. All the potentials on the 

inverse boundary can then be calculated and, since the span- 
wise and vertical derivatives are small, will primarily be func- 
tions of the pressure coefficient at grid point i — X A and the 
value of the potential at grid point /— 1. 



Surface Calculations 

As the inverse boundary conditions drive the flowfield to a 
converged solution, it is necessary to calculate periodically the 
location of the new displacement surface and to regenerate the 
computational grid about this new geometry so that the 
pressure boundary surface will correspond to the physical 
boundary surface being designed. Each new surface is found 
relative to the previous surface from an integration of the wing 
surface slopes. The surface slopes are calculated from the cur- 
rent flowfield solution using the flow tangency boundary con- 
dition, which in curvilinear coordinates is 

V 1 x VS = 0 (15) 

Note this condition, with the gradient in the curvilinear plane, 
is a direct analog to the same condition expressed in the 
physical plane. 

A more useful expression can be obtained by expanding the 
above equation to: 


( dv > 

) — 

- — 1 

fj2_) 

\ d( J 

U 

u 1 

l aj- / 


This expression can be solved for the new chordwise airfoil 
slope <h)/d£ if the current values of the spanwise slope 
are used. Since the wing surface is represented in the computa- 
tional grid as a plane of constant t?, the current slopes on the 
wing surface equal zero and a simplified flow tangency condi- 
tion results: 


(dn/a{) wini = F/t/ (17) 

The above expession has been applied to the computational 
surface plane in order to find the relative location of the new 
physical surface. This approach is an approximation, since the 
above equation is only exactly true when applied to the new 
surface itself. Using this method, however, provides for a 
stable iterative surface updating procedure that quickly con- 
verges 6 to the target surface. 

To calculate the relative surface slopes, it is first necessary 
to determine accurately the values of the contravariant 
velocities, U and V . As was also determined by the work of 
Weed et ah, 12 a simple finite-difference calculation of these 
velocities is not sufficiently accurate. Borrowing from Ref. 12, 
a more accurate method was implemented that uses the 
residual expression to calculate the velocity ratio V/U , under 
the assumption that the residual is zero at the surface points. 
The residual expression from FLO30 can be written in finite- 
volume form as 

+ (other terms) * 0 (18) 

The“other terms” in Eq. (18) involve the grid point coupl- 
ing and upwind dependence terms of the formulation and are 
assumed to be constants in the following development. 

The desired velocities can also be written in this finite 
volume form as: 

V phV f^phV) 

V fihU n^iphU) 1 y) 

By simple manipulations, the normal velocity can be ob- 
tained from the residual expression as 

{phV) =2fi n (p/iK) f _, (phU) 

— (ph W) — (other terms) (20) 

where the subscript ij - 1 refers to the values at grid cell centers 
above the wing surface. 
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automates the iterative selection of a leading-edge shape 
would otherwise have to be performed by the user. 

Two method* of relofting can presently be seized. The 
first method is a simple linear rotation scheme. This m 
can be visualized with the help of Fig. 6. The dashed line m- 
dicates the original leading-edge geometry and a hypotheucal 
new surface shape that has been calculated for the mverse 
design regions. Without modification, this new surface has a 
traihng-edge thickness of A. If a thickness of A » " ere s ^ r, * d 
by the user, then the surface would have to be relofted o 
changed. In the present scheme, in order to obtain the desired 
thickness, a displacement thickness 6, is i added to ihecurrent 
design surface. This thickness has a distribution from the 
leading to the trailing-edge and is determined by the formu a 


i.(x) = (A,— A)(x/c) 


( 21 ) 


Fig . 5 Comparison of slope calculaUon methods. 

In order to use Eq. (20) to find the desired surface velocity 
ratio, it is necessary to know the U and W velocity com- 
ponents at the "ghost” cell centers below the wing surface. 
These values can be obtained in a manner consistent with 
FLO30 by specifying the “ghost” cell values to equal Uie 
values at corresponding points immediately above the wing 
surface. A comparison of the accuracy of both the finite- 
difference approach and residual approach is shown in Fig. 5. 
The calculated displacements are for a converged analysis 
solution for which the calculated slopes should, of course, be 

With the contravariant velocities known, an integration of 
Ea. (17) in the chordwise direction £ from the start of the in- 
verse region to the trailing edge yields a set of surface 
displacements for the new wing surface relative to the previous 
one. These displacements are expressed as changes in the com- 
putational coordinate v and are converted to surface 
displacements in the physical plane 6{x) via the local gnd 
transformation. The physical plane displacements are coinci- 
dent with the computational grid points in the inverse regions. 
To obtain the corresponding displacements at the original 
geometrical locations specified in the program input data, a 
Unear interpolation of the above data is performed. Finding 
the displacements at the original geometry stations permits the 
calculation of the new wing airfoil sections at the same 
semispan locations. 

Trmiling-Edge Closure 

The procedures outlined above will compute a wing surface 
corresponding to a given, fixed, leading-edge geometry and to 
a desired set of pressure distributions in the inverse regions. 
The above procedures do not, however, guarantee that this 
wing geometry will be practical. In particular, past experience 
has shown that inverse surface calculations may yield airfoil 
sections that have either excessively blunt trailing edges or 
which, at least numerically, have the upper and lower surface 
crossed at the traiUng edge (“fish tailed”). The former case is 
undesirable due to aerodynamic considerations, while the lat- 
ter is physically impossible and may produce unpredictable 
problems in the grid generation or flow calculauon portions of 

FI since for any specified pressure distribution the correspond- 
ing wing surface will be controlled by the leading-edge 
geometry, which serves as an initial spatial boundary condi- 
tion for the inverse region, the problem of assuring trailing- 
edge closure can be viewed as the proper selection of a leading- 
edge shape. A procedure for systematically modifying the 
leading-edge region in order to achieve some desired traihng- 
edge thickness is called relofting. Such a relofting procedure 
has been incorporated into the present design process in order 
both to prevent the problems of trailing-edge crossover and to 
allow the user the option of specifying a trailing-edge 
thickness as an additional design variable. This design feature 
should be very useful in practical applications since it 


where c is the chord length of the local airfoil section. The 
total displacement for a surface update is then the sum of the 
two displacements 6(x) and «.(*). When both the u PI* r “ d 
lower surfaces are designed simultaneously, the displacement 
magnitudes determined by relofting are divided between the 
two surfaces so that half is added to the lower surface and half 

to the upper surface. , „ 

The second relofting method uses the same approach as the 
first for the aft inverse regions, but modifies the 
region by a proportional thining or thickening of the surface 
ordinates. This approach can be expressed by: 


y* x (X) =y",*'\y n {x)/y”j] 


( 22 ) 


here they subscript refers to the ordinate at the direct-inverse 
junction determined from the linear relofting ; of the ^aft 
regions. Note that this method will produce leading edges in 
[hesame family of shapes and, for example, allow the design 
of an NACA 0006 airfoil when starting from an NACA OUiz 
airfoil (sec test case II). 

Results 

A variety of different test cases were run as verifiration of 
the current design method. These cases involved both sub- 
critical design and supercritical design over secuon 6 e0 ‘™ :t " 
selected to test the versatility of the input and design comro 
logic. In this section, results from two of the more significant 
test cases will be presented. The results shown wcre °btamed 
on a medium grid having 81 streamwise, ^ vertical, and 19 
spanwise points with 11 spanwise stations and 53 points on the 
wing at each station; and, in all cases, the maximum change in 
the reduced potential was reduced at least three or d= f s 
magnitude. Thus, the results do not represent ultimate <»n- 
vergence but should be representauve of engineering 

accuracy.” , . 

The use of the medium grid for the design cases shown m 
the following was dictated by computational cost and tune. 
Fine grid solutions for these type geometries have been ob- 
tained but are not significantly different from the medium pid 
results except for a generally smoother shape. Use of the fine 
grid in design is necessary, however, when the airfoil «rtions 
involved are aft cambered, since a higher gnd-point resoluuon 
is needed in the trailing-edge regions. . 

The planform selected for the test cases was the Lockheed 
wing A wing-body. The wing for this configuration has a 
quarterchord sweep of 25 deg, a Unear twist distribution rang- 
ing from 2.28 deg at the wing body junction to 2.04 deg at 


Original Design Surface 



Fig. 6 Trailing-edge thickness ndjusted by relofting. 
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the wing tip, an aspect ratio of 8, and a taper ratio of 0.4. The 
last two values are based on the wing without fuselage. 
However, instead of the supercritical sections normally 
associated with wing A, the initial airfoil sections at each span 
station were assumed to be composed of symmetric NACA 
four-digit airfoil sections. 

The target pressure distributions used in the design regions 
for the first test case were selected to yield airfoil shapes 
thicker in the aft portions of each section and, at supercritical 
conditions, to yield on the upper surface weaker and more for- 
ward shock waves than those that would normally occur on an 
NACA 0012 section. On the lower surface, the target pressure 
distributions were selected to have either a favorable pressure 
gradient or fairly constant pressure plateau over much of the 
lower surface. 

For the second test case, the pressure distribution was ob- 
tained from analysis solutions of an assumed wing geometry. 
The intent of this case is to verify the relofting procedures and 
show the ability of the current method to make large surface 
changes in going from a thick wing to a thin wing (approx- 
imately 12% to 6% thick, respectively). 

Both cases were for a freestream Mach number of 0.8 and 
an angle of attack of 2 deg. In each case, the pressure distribu- 
tion was specified in the design regions from the 15% local 
chord location to the trailing edge and used as the boundary 
condition in these inverse regions starting with the first itera- 
tion. Prior to the first design surface update calculation, 300 
SLOR iterations were executed and, subsequently, surface up- 
dates were computed every 50 cycles. The solution was con- 
sidered converged and terminated after 550 total iterations for 
the first case and, due to the large amount of relofting re- 
quired, after 950 iterations for the second case. 

Test Case I 

The inverse design regions for case I, which was an attempt 
to design both upper and lower surfaces on two noncon- 
tiguous regions of the wing at supercritical conditions, are 
shown in Fig. 7. A comparison between the initial pressure 
distribution associated with NACA 0012 sections and the 
target pressures for two of the designed sections is portrayed 
in Figs. 8 and 9. As can be seen, the target pressure distribu- 
tion essentially eliminates the upper-surface shock wave pre- 
sent at inboard stations of the original wing; at outboard sta- 
tions, it weakens the shock and moves it forward. In addition, 
significant changes in the lower-surface pressure gradients are 
evident. Also shown in Figs. 10 and 11 are the pressures com- 
puted by the program at the end of the inverse design pro- 
cedure (denoted as M design pressures”). These pressures are in 
excellent agreement with the target pressures, which indicates 
that the method is satisfying properly the desired inverse 
boundary conditions. 

The corresponding designed airfoil sections for this case are 
shown in Figs. 10 and 11. Even on the expanded scale, the 
agreement between the designed and target surfaces is ex- 
cellent at all design stations. However, trailing-edge closure 




(cmsc 1, 30% station). 



(case I, 70% station). 



was not enforced for this case, and at the boundary stations 
there is some departure between the designed surfaces and the 
target surfaces near the trailing edge. It is believed that this 
slight difference is a ramification of the change in spanwise 
regions. 

In any event, the pressure distributions resulting from an 
analysis of the designed surface shown in Figs. 10 and 1 1 are in 
excellent agreement with the target pressures, as can be seen in 
Figs. 12 and 13. In addition, the section lift coefficients at the 
various design stations are in very good agreement with the 
target coefficients. Based on these results, it is believed that 
the present method can adequately design/modify nonadja- 
cent regions of a wing in transonic flow. 
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sections (ease I, 70% station). 



Test Case II 

This test case was selected to demonstrate the ability of the 
design methodology to handle two difficult design tasks. The 
first task was to change a wing from supercritical to sub- 
critical, which is both a typical engineering task and a signifi- 
cant problem for wing design algorithms. The second task was 
to make large surface changes to the original airfoil without 
generating large surface distortions from the accumulation of 
geometry calculation errors. Due to the upstream dependence 
of the supersonic flow, this required making large changes in 
the leading-edge region through the rclofting procedures. The 
design regions for this case are shown in Fig. 14, where the 
wing thickness varied from 12 to 6% between the wing root 
and 80% span location, and was constant going outward to 
the tip. The input design pressures were for a constant 6% 
thick wing. 

The first attempts at this design used the linear leading-edge 
relofting procedure and from a practical standpoint were un- 
successful. The final design surfaces were still supersonic in 
the leading-edge regions while satisfying the subsonic aft sur- 
face conditions by producing strong shocks at the direct- 
inverse junction location. In addition, the surfaces themselves 
had sharp surface slope discontinuities at the same location. 

When the thinning approach was used to reloft the leading 
edge, much better solutions were obtained. Figures 15-26 
show the changes in pressure distribution and surface shapes 
with a comparison of target to designed surface pressures for a 
few span sections as in the previous case. As can be seen, ex- 
cellent agreement between target and final pressures and sur- 
face were again attained for this extreme case. The only 
noticeable surface irregularity are a small wiggle at the direct- 
inverse junction that can be seen as a small pressure wiggle in 
the pressure plots. 



Fig. 13 Comparison of pressures from analysis of designed wing 


with target pressures (case I, 70% station). 



Fig. 14 Design case II. 




Fig. 16 Comparison of initial pressures with target and final values 
(case II, 20% station). 
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Fig. 17 Comparison of initial pressures with target and final values 
(case II, 40 ft station). 



(case II, 60ft station). 



Fig. 19 Comparison of designed sections with original and target 
sections (case M, Oft station). 




Fig. 21 Comparison of designed sections with original and target 
sections (case 11, 40ft station). 



sections (case 11, 60ft station). 



Fig. 23 Comparison of pressures from analysis of designed wing 
with target pressures (case Q, Oft station). 



Fig. 24 Comparison of pressures from analysis of designed wing 
with target pressures (case 11, 20ft station). 






NOVEMBER 1988 


TRANSONIC WING DESIGN 


1017 



Fig. 25 Comparison of pressures from analysis of designed wing 
with target pressures (case II, station). 


- 1.0 * 


-0.5 H 


o o.o H 



design case n 

60* SEMI-SPAN STATION 

TARGET Ci - 0.153 
DESIGN C i -0.151 


DESIGNED SURFACE 
PRESSURES 
TARGET PRESSURES 



Fig. 26 Comparison of pressures from analysis of designed wing 
with target pressures (case II, 6097* station). 


Conclusions 

A direct-inverse wing design method has been successfully 
incorporated into the TAWFIVE transonic wing-body 
analysis computer code. The resultant code is capable of 
designing or modifying wings at both transonic and subsonic 
conditions and includes the effects of wing-body interactions. 
A series of test cases have been presented that demonstrate the 
accuracy and versatility of this inverse method. Additional test 
cases and results axe also presented in Refs. 16 and 17. 
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Abstract 

Progress in the direct-inverse wing design method in curvilinear 
coordinate* ha* been made. A «par.wi.e oteiUaOon problem and pro- 
posed remedies are discussed. Test eases are presented whrch reveal the 
approximate limits on the wing’s aspect ratio and leading edge wing 
sweep angle or a successful design, and which show the s.gnjficance of 
span wise god skewness, grid refinement, viscous interaction, the initial 
airfoil section and Mach number - pressure distribution compatibility 
on the final design. Furthermore, preliminary results are shown which 
indicate that it is feasible to successfully design a region of the wing 
which begins aft of the leading edge and lerminates prior to the trailing 
edge. 


Introduction 

With the advent of efficient numerical schemes that accurately 
model the irrotational transonic flow about complex configuration* *uch 
as wing-bodies and the appearance of computers with memory capac- 
ities and computational speeds necessary to execute these schemes in 
a reasonable amount of time, the efficient design of wings for tran- 
sonic flight is quickly becoming a reality. Although transonic potential 
schemes combined with integral boundary layer solvers may not model 
the real flowfield as accurately as Euler or Navier Stokes Schemes, their 
use can significantly reduce the costs and time expenditures associated 
with transonic wing design. 

Many methods ranging from optimisation techniques 1 - 2 to var- 
ious inverse methods have been formulated using potential solvers to 
design wings in transonic flight 3-12 . One such method, which has been 
under development at Texas A£:M University for the last several years, 
is the direct-inverse transonic wing design method. In this method the 
airfoil sections making up the wing are created by specifying desired 
pressure distributions over all or part of the wing aft of the leading edge, 
solving via finite-difference or finite-volume techniques the mixed Neu- 
mann and Dirichlet boundary value problem associated with the full 
potential equation for compressible flow, and then integrating the flow 
boundary condition at each spanwise station in the design region. The 
design pressure distributions can be selected by an experienced user 
to have such desirable characteristics as mild or nonexistent shocks, a 
slowlj r increasing adverse pressure gradient, a center of pressure giving 
a desirable pitching moment, or an efficient spanwise loading, The de- 
signer may also use wind tunnel tests of successful airfoils as an aid in 
picking .a desirable pressure distribution. 

The direct-inverse technique has been successfully used in a 
stretched and sheared cartesian system*"* and more recently in a 
curvilinear system. 9 ' 11 This paper presents progress in the latter. It' 
will include a brief description of the analysis and design methods, 
techniques used to suppress a spanwise oscillation problem resulting 
from the interaction of the design method with the potential solver, and 
a scries of test cases. The latter will reveal the lack of dependency of the 
wing design on the initial airfoil section, the importance of including 
viscous effects in wing design, and constraints due to aspect ratio, wing 
sweep and spanwise grid skewness. 

Background 

FLQ30 

The original computer program, which was modified intinlly by 
Gully 9-11 for inverse wing design, is TAWFIVE 1 - 13 . This program* 
not only has the capability of computing the potential flowfield about 
c. fairly general wing and fuselage combination, but also contains a 
three dimensional integral boundary layer scheme which provides the 
necessary viscous effects in the form of the boundary layer displacement 
thickness, wake curvature and wake thickness 13 . 


The inviscid numerical scheme is based upon Jameson and 
Caughey’s conservative, finite-volume, full potential flow solver, 
FLO30, in which computations are performed on a body-fitted, sheared, 
parabolic, wind-tunnel type coordinate system. The theory behind this 
code will only be briefly discussed here, but further details can be found 
in Refs. (14-J8). 

FLO30 solves the compressible potential flow equation in conser- 
vative form written in curvilinear coordinates 

lf>hU) ( +(phV) v + (ph\V) ( = 0 (1) 

where the nondimensionalised physical velocities, u, v, u/ are related to 
the gradient of the reduced potential function, 4 > , hy 



and the nondimensionalired eontravariant velocities U t V, W are related 
to the physical velocities by 



Here H is the transformation matrix defined as 


*< 

V( V*> VC 
H £ ( 


with h = |H | 


( 4 ) 


The local air density normalised by the freestream density can be 
obtained in terms of the local speed of sound from the energy equation 
and isentropic relations as 


p = (sKc,)^ 5 


( 5 ) 


Hd the local speed of sound normalised by the freestream 

can be calculated in terms of the nondimensionabsed speed of sound a 

stagnation, a«, using 


a==a!-f-2-Vfe^ 


(61 


Flo30 uses a finite-volume type scheme which makes use of a staggered 
box approach. Us formulation is directly analagous to the control 
volume approach used to derive the original PDE in Eq. (1), except 
in the finite-volume scheme, the discrete nature of the finite difference 
model is considered from the onset by using a finite control volume in 
the neighborhood of a grid point in the finite-difference mesh 19 . This 
method is best explained by using it to discrete the following two- 
dimensional, incompressible version of Eq. (1) written in cartesian 
coordinates: 

= 0 ( 7 ) 

with the aid of the two-dimensional box shown in Fig. 1. As can be seen, 
from the figure, the staggered box scheme derives iU name from the 
way in which the primary and secondary boxes interlock. The values of; 
the potentials at the four grid points which make up the corners of each; 
primary box are used to calculate the velocities, u, r, in the following 
manner: ' ' 

u = 4> x =Vy6 t 4> ^ 

v = 4>y =Pz&y4> 
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where p and S are averaging and differentiating operators respectively 
and are defined by Jameson as: 



M*/ = o (9) 

— /i+jj — 

• • J ,K«t A* = 1. Therefore, the velocity, u, for 

where it is assumed that A* i y , • found bv 

instance, at the primary box center located at (« + ,.J + 3 J 

'*, + !.j+l + l (10) 

f<* , -d. .1+ zAzlli 

2 

, k ..,-r.ndarv box is determined by 
The flux at the midpoint «' h ^ lhe £orners of that box in the 

averaging the velocities, u ' „ inl0 tb e secondary box 

. Ltfl-sssrssSiiS-' 

p,i, W + Mv M = ° (U) 

where for example 

(u. + 1 i 




( 12 ) 


.tatement of the conversation of mass and should tend to zero as the 

""tSS earlier, a body-fitted, wind t»n„«l-typ« grid is used 
in FL030. The grid shown in Fig. 2 is the coarsest meshandha.40 x 6 
x 8 point. in the {. V. *"d < direction, respectively. V, . »th this the 

wing becomes a constant ij surface, and each cylindrical looking 
i, a constant < surface. Constant < lines can be seen running span* r se 
on the wing at constant chord fractions from the leading edge Not.ee 
also due to the conformal transformation used , that constant f lin 
are ’packed close to the leading edge of the wing. This ‘lustenng , 
an attractive feature when designing airfoil sections using the direct- 
inverse approach. Moreover, constant < lines are spaced evenly on the 
wing and, on the finest mesh, give the designer up to -1 sp.nw.se 
stations where the pressure distributions can be specified. 

The flow tangency boundary conditions are enforced by reflecting 
the fll« abovetm surface to ghos, points below the wing, use «£. o 
symmetry plane such that the net out of plane component of the flux 
iero at the surface. 
jpverse Design Method 

In the direct-inverse method a pressure boundary condition is 
enforced in the design region rather than flow tangency. As *hown b> 
Gaily’’' 10 , the input pressure coefficient can be written in term, of the 


m the previous discussion, the OTplicitjassumption 

the°flux hnto'the 11 top'of thetseconday cel. face wou.d be, for instance: 


fr. ,t-i m if 

r-jL)] 

-r — 1 

- 1 

[[i + “ M “ V 

9 = oo/. 



(19) 


/ * v{x,y)dz~ 






+ ‘ 

T 


(13) 


where q 2 = (v 2 + v 2 + ti> 3 ) (Too • 

Solving for « in Eq. (2) and Eq. (19) and then equating the two 
results give* 

Jn 4>t + = 




_ , Cnuehev 15 found that lumping the fluxes at the primary 

Jameson and Caughey solution between adjacent grid 

cell centers led to an uncoupling of the ^oluuo basically 

tzz «—»-«*« ‘ fi,i 1 ““ i “ 

(i-rl.J + i) wou,d be 

25 f—) (M) 

■i+ij+j Vey/.»j+i 
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(20) 


where Jij are the elements of (h T ) 

A potential, can be found in terms of the pressure coefficient 

v 1 • p_ at the grid point location (t ;>•?» )» an 

by applying Eq. (20) at gn P . and seeo „d order 

using central differences in the (. and < direction an 

backward differences in the normal direction, ij, yie g 




where » 

( — 1 W+l 

(15) 

1 

\d Vi+lj+lJ 


+^r,kj-2,t "h j 

When all the fluxes are extrapolated in this manner 

and included in 

+^r-l,ky,k + l 
4 


Eq. (11), the result is 

/I,,*..* + H** 1 ,!* ~ ~ 

, aDDlied to a three dimensional flowfield in 

> T«“ (A® . P) * « A IAV 9) • + «.<< t*w + « 

/ r /I ) — 0 

where the Q’s are the compensation terms defined as 

£?<, = (A ( + A,)rW 

Q„< = (A, + A< )Mt«u< * ( 18 ) 

<3 a = (A< + A ( )jiA 

Q (£ = (A { +A„ + A ( )d(e<^ 

. . . - te i n nuence coefficients which compensate for the depen- 

A ( ,A„A ( is, is p Q a nd H are upwind.ng terms 

dence of p on an exclude un- 

which desymmetrire the sche, " e ^ S " P "*° tion by providing an arti- 
real, discontinous expansions fror " tht ' OR f s of cour se a direct 
ficial viscosity. This equation, solved via 5LOK, 


+ E 


{ C > 


( 21 ) 


where F (c, b the right hand side of Eq. (20) and y = h y on 

th ‘ ilnce’thlTrid is boundary conforming, the wing £ the 

desigm region must be updated every so *-££*** ^ ^ 
tangency condition written in curvilinear coordinates 

( 22 ) 


(dr,\ _Y. - — 


The integration of this equation can be handled in two different ways 
K the spanwise term, g. is lagged one g.oba! iteration, it wil alwa 

respect to the ( or < direction vanish on the wing's surface, and, Eq. (22) 
reduces to „ . x /t/ \ 


fM = (1) 

\^W (.»»,» \U / i,ky,k 
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Tht other approach would be lo integrate Eq. (22) iteratively. If 
the contravariant velocities are froien at their current values, and 
the spanwise terms are initially assumed to be sero, Eq. (23) can be 
integrated lo find the approximate inverse changes At). These can 
then be used to find approximations to the spatial spanwise derivative, 
These can then be included in Eq. (22) to provide a better 
approximation to the flow tangency equation, and the process can 
be repeated using Eq. (22) until the spatial derivatives are converged. 
Numerical experiments reveal that the spanwise terms are at least two 
orders of magnitude smaller than the chordwise terms pnor to the 
creation of the new grid. Hence, the spanwise terms can be normally 

n eg IC T*ht contravariant velocity, V, can be obtained most accurately 
from the residual expression* 10 . Combining the previously defined; 
averaging and differencing operator* 

M, (pM'),s,.s = i (OU'Wm + »•■) (24) • 

«, (phn.s,.s = {(p hV kk,-\.k - (!5) 


yields 


i, (phV) t ^ k = 2 - 2^ (pb^),.s,.i 


(26) 


Substituting this into the residual expression, Eq. (1<), and solving for 
the out of plane flux, phV, on the wing surface yields 


2pt>i( (phV);,i,,i = Pti<G( (phP)i,s,,s + 2p(( 

+ MtA (pW.,*»,») 

-f compensation and upwinding terms 



Since at convergence the flow should be tangent to the designed surface, 
the tangency condition is enforced in the residual expression, Eq. (10. 


by setting 

(pM')i. !,,+}.* 


(28) 


Since the resulting expression is identical to the RHS of Eq. (2i), the 
expression for the normal fiux becomes 


Residual 

(P^Ot.ky.k - ^ 


(29) 


Note that since the residual obtained in the iterative process is not 
initially zero in the design region due to the application of the inverse 
boundary condition, Eq. (29) reveals that there will be an ejection of 
fluid from the wing boundary. 3 - 4 ' 3 * No attempt was made in the present 
iterative scheme to account for this temporary addition of mass into the 
flowfield, since at convergence it would be negligible. Upon substitution 
of Eq. (29) into Eq. (23) and using the cell averaged fiux phU on the 
surface, the equation for the required slope change becomes 


djn _ _ Pin((P kV ) _ R' sidual 

di ” V % p iv< {pW) 2 (pW) 


(30) 


The changes normal to the surface at each spanwise station are obtained 
by integrating from the beginning of the inverse region to the trailing 
edge using the trapesoidal rule. (Higher order integration schemes were 
tried but had little effect on the Anal answers, except for coarse grids 
in regions of high curvature such as the cove region of a supercritical 
airfoil.) Assuming that the grid line leaving the wing in the t? direction 
is normal to the wing, these changes, Ai), are then converted from 
computational lo physical units by scaling them by the transformation 
metrics such that 


A l =s At) 



+ 


ay’ 

dr) 


(31) 


After subtracting the boundary layer displacement thicknesses from 
the inverse corrections, At, which have been linearly interpolated to 
the user defined input stations, the resulting changes are added to the 


initial airfoil sections yielding the new wing surface for the current time 

'' Vtl Many times, the trailing edge thickness may be too large if the 
leading edge curvature is too small or may be ‘fish-tailed if the leading 
edge curvature is too Urge. These undesirable situations are remedied 
by a procedure called relofting where the designed surface is rotated 
about the leading edge to meet a specified trailing edge ordinate. • • 

This relofting procedure*-" is usually carried out for every surface 
update. To illustrate the previous procedures, the first global iteration 
of a typical design before and after relofting is shown in Fig. 3. 

Spanwise O scillations 

An annoving divergent spanwise oscillation problem sometimes 
occurred when designing a wing which required extensive rdoft.ng 
especially when the initial section was thinner than the target This 
oscillation led to sections which were too thick or too thin at adjacent 
constant C 8^ *ftionf. An example of this phenomenon is shown 
in Fig. <■ This problem was more pronounced when the sweep was 
increased or the aspect ratio was decreased. After many failed attempts 
at remedying this problem by reformulating the inverse boundary 
condition, attention was directed towards the residual and the term 
composing it. The residual was broken into its components and plotted 
after each surface update. Sample plots for a divergent ease are shown 
at four different time levels in Fig 5. As can be seen, the compensation 
terms that include spanwise derivatives of 4 are at first very smaU 
compared to the rest of the terms, but later tend to dominate and 
amplify the oscillation. This oscillation starts at the direct-inverse 
interface or. in other words, at the first spanwise station from the root 
in the design region mid propagates spanwise as adamped o«ullaUon 
with a period of two grid spacings. Presently, it is believed that the 
initial mismatch in the potentials at the direct-inverse tnterfaee in the 
spanwise direction is amplified by those compensation terms which 
include spanwise derivatives of the potential function. The res.dual 
is then undershot and overshot on alternating spanwise stations. This 
oscillation is further magnified by relofting, which creates a section 
that is too thin when the slopes defined in Eq. (30), which of course nre 
directly proportional to the residual, are too largeandvce-versaSince 
more or less fluid has to be ejected from the section that is too thin 
thick respectively to give the streamline approximately correspond ng 

£1 SS ui... La-, a. r— au «•“ •* •“ l d “‘” 

i, taken further away from the adjacent fields by the inverse boundary 
condition which forces an even further unoershoot or overshoot of the 
residual. It should be noted that this problem is not soley due to the 
implementation of the direct-inverse technique since this oscillation as 
not been observed with the ZEBRA design code’ but seems l ° *' 
unique to the coupling of the method with the analysis code FLO30. 

After exploring many alternatives lo counter this problem, four 
methods have been developed to damp out the spanwise oscillation: 

A) Specify the inverse boundary condition at at least every other 
spanwise station and linearly interpolate the inverse displacements to 
the stations lying in between. This has been named the type 11-2 

method. ^ j nverse boundary condition at every station, but 

again only calculate inverse changes at every other station and linearly 
interpolate the inverse changes to the stations in between. This will be 
referred to as the type 11 method. 

C) Immediately prior to every surface update, calculate all span- 
wise derivatives of the potential used in the residual based upon a poten- 
tial function smoothed in the spanwise direction. This is accomplished 
by first defining tbe smoothing operator as 

0 < £ < 1 (32) 


<r c / = “ j) A + 4 A+i» 


where c determines the amount of smoothing. Then using c in the 
spanwise differentiation of * with the maximum amount of smoothing 

(i.e., t = 1) 


f—) 




the smoothed spanwise derivative of 4> becomes 


. _ (jjjA+i - + ti .A** 1 

(♦cJijA+i “ 4.0 


di,jk-i) 


(33) 


( 34 ) 
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oscillation but only slowed At ^ ^ thf fout different approaches. 

Tnd a 1 Mach 

iniual sect , op »n the *^ re ^ &djng ed ge B nd extended to the 
semispan and began 5 /o aDDT0Be hes worked well for this case, by 

trailing edge. Although all o PP . 7 between the target and the 

using the A and C produced 

designed sections as a me J we „ „ Bl lhe edges of the design 

the best results in the , ^ nuI „ber of flowfield iterations, 

region. On the other ha , unsal ; jfBCtory results when compared 
method D produced approach on the residual at 

to the target sections The £ of t PP ^ ^ p . g g 

the trailing edge after 10 surl P . due l0 the fact that 

discontinuities in the residual fot tnelho, J0 50 an d 70% 

,h, i.... mm J»r » 

in the compensation terms ^ the fouI methods (not 

, rr^s^hat the smoothing approaches (methods 
r^d D1 work well when designing in the interior of the wing, but di 
„°t give large errors 

typt II ..a W '« -W. -* 

- ttflSSi =ss 'z^z^;:£z 

S£Sr^s£«« 

the quickest. 

Results 

a~ u. «»*>;: £ ilrifs,— .n; 

inverse design procedure. The case we e ^ ^ ^ ftn(J 

imate limits viscon* interaction, grid refinement, 

leading edge and terminating forward of the trailing edge. 


.Span wise Grid Skewness 

Recently it was discovered that the skewness of the constant ; ( 
grid lines leaving the tip of the wing (Fig. 9 ) can have a dramatic 
efTect on the design of the sections near the wing tip As ““ *« 

in Fig. 10. if the design grid was significantly skewed and the input 
pressures were obtained from an analysis on on unskewed gnd, the 
converged design yielded incorrect airfoil shapes in the tip region. This 
difficulty is due to the large errors near the wing tip associated with a 
skewed grid which are revealed in the pressure distributions ( lg- )• 
The grid skewness has caused the shock location to move further aft 
Although the skewness of the grid was quite extreme in this case, these 


possible. 

Boundary Laver and Wake pffecU 

. . .. f il;. *tudv was to determine tbe sig 

One of the objectives '*£ of transonic wings-. 

nificance of various viscous tKe ^ "^J hvtion WBS obtained by 
To study these effects, an input pr«s included 

analysing and wake cur- 

boundary layer displacement thickness, degrees and 

rri7:‘iS‘^ =a. » - *-»? * • «■- 

computations were performed on a fine ( 160s ^?-J ln the 

This pressure distribution wi^Jien ased h ^ the 

first case, the w.ng was d«.gned .nv,sc,d boundary 

wing was designed without the wake option but .n«^ ^ ^ 

displacement thickness elr * ct *- " 1 design region for all three 

were used in the design of the wing. T . h ' d “' 8 " of lhe lading 

-a 

edge of the “ ri “ L . the 30 5 0 and 70% semispan stations; 

=^s£S553=SS 

were obtained through linear interpolation Neglecting 

1 

for this, the target wing was analyx reVeB l that at 

The resulting pressure distnbufons show « t«F^ ^ & very sm all 
the Re chosen, wake curvature an rt distributions, 

- ,b, Mi -a W. it ta*, 

except near the trailing e g ; n vestieated it discovered 

layer displacement thicknesses we ^ ^ p rod u C ed boundary layer 

SSSSoit-' « -• » Ji% 

initial inverse changes to yield the hard airfoil, these larger ^placement 
thicknesses would produce a section that was thinner than the target, 
but, after relofling the airfoil section would actually be thicker than 

thC ‘“he wing sections designed inviscidly are profoundly- dlffemnUt ‘ 
differences at the inboard and outboard design stations are due to th 

clinwn in p: c i 4 led to airfoil sections which vaned smoothly in 
case, shown in g. , - t . t ; ons In Fig. 14, the boundary layer 

the spanwise direction at all stations, in r g , lareet 

displacement thicknesses obtained from a viscous 

were added to the target sections for comparison with the invisodiy 

designed sections. , 

' After the wings were designed, all three were then analysed w it h 
full viscous effects to assess the significance of the changes made to the 
wing on the pressure distributions and to see how well these P'“ su f** 
matched the target pressures. Knowing that the wmg im** wj 
full viscous effects is correct, it is quite obvious from Fig. 15 and Table 
1 that the wing designed inviscidly is quite unsatisfactory. The shock 
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is not far enough aft and the lift produced is sometimes 20% smaller 

than g^ t e j” , " lhe rc4u | ls of this study, it can be concluded that for 
the Re and Mach number chosen, w.ke curvature and wake th.ckn«s 
have a very small effect on the airfoil sections designed. However, the 
boundary lavei displacement efTect has a profound effect on the section 
shapes and hence must be included in the design process to yield a 
wing which will produce the desired lift in a viscous environment. 

Wing PJa pform Effects 

Three cases were attempted to roughly delimit the applicable 
range of aspect ratios and leading edge sweep angles. These included 
Lockheed Wing A, B and C. These wings have aspect ratios of 8, 3.8, 
and 2.6, leading edge sweep angles of 27, 35. and 4» degrees and taper 
ratios of .4, A, and .3 respectively. The target pressure distributions 
were obtained by a direct analysis of the target wings in an inv.scid 
environment. The initial section for Wing-A was a NACA 001., while 
a NA.C4 0006 was used for Wing-B. The original section was used with 
Wing-C due to the difficulty of the case. Also for Wing-C, as opposed 
to the circular cross-section, an elliptical cross section of the fuselage 
was used to provide a flatter surface for the grid generation package 
The circular cross-section combined with the large relative t ic ness o 
the root section compared with the width of the fuselage played havoc 
on the grid at the root, as can be seen in Fig. 16. 

In order to better understand flow about each wing, the corre- 
sponding velocity vectors on the surface of each wing were plotted, as 
shown in Figs. 17-18. As to be expected, the spanwise component of 
the flow increases as the aspect ratio decreases and sweep Incre “«- 11 
is also interesting that there seems to be an inboard flow for all three 
cases on the upper surfaces aft of the leading edge. If this is correct 
this inboard flow may be attributed to the effect of the fuselage and 
the wing tip vortex. 

The design region for Wing-A wd Wing-B extended from 10- 
100% semispan and began 5% and 2.5% aft of the fading ed 6 e ' 
respectively. Computations were performed on a fine grid. Results 
for Wing-A are shown in Fig. 19, while results for Wing-B are shown: 
in Fig. 20. As can be seen the designed and target sections for both; 
wings are in excellent agreement in the interior of the design region and 
closely match at the edges of the design region. 

In the case of Wing-C, the section shapes should not have changed 
with the application of the inverse boundary condition. But, because of 
the large amount of spanwise flow and the associated spanwise gradients 
for Wing-C, the spanwise oscillation effect could not be overcome with 
any of the present remedies. Further information about this case w as 
obtained by using the Type II method and not relofting the section 
shapes. The results for such a converging fine grid case are shown in 
Fig. 21. The first design station at 18% semispan is too thick on the 
upper surface as compared to the target. This discrepancy is again 
due to the over prediction of the residual at the first station due to the 
initial mismatch in the potentials in the spanwise direction, and, hence, 
to large spanwise gradients of the potential. The errors dimmish as the 
tip is approached, but are always relatively large in the trailing edge 
region due to the difficulty in accurately imposing the inverse boundary 
condition near the trailing edge for this case. If an attempt were made 
to converge this case further by continuously relofting the shapes to- 
meet the trailing edge ordinate, the same spanwise oscillation problem 
would again occur. 


Initial Profile Effects 

One of the disadvantages of the direct-inverse method is that a- 
priori knowledge about the shape of the leading edge must be known 
to achieve suitable airfoil shapes and desired trailing edge thicknesses. 
Relofting does alleviate this to a large degree; but it will not, in general, 
happen to produce a leading edge that will yield the desired pressure 
distribution in the leading edge region if the inverse boundary condition 
is by necessity applied too far aft. It was thought that because FLO30 s 
grid package dusters grid lines close to the leading edge of the airfoil, 
that the design could be started quite dose to the leading edge, thus 
relieving the designer of the difficulty of choosing a correct nose shape. 
Two test cases were conducted to investigate the dependence of the 
final design on the initial airfoil section. Both used Lockheed Wing-A 
at the same conditions mentioned earlier for the viscous study. For t e 
first ease, the initial airfoils were the same as those in the viscous study. 
These airfoils all had leading edges which were in the same family as the 
target section. The design was started 10% aft of the leading edge. In 


the second case, NACA 0012 sections were used at all design stations; 
the leading edge of these sections were not in the same family as the 
target airfoil sections. For this case, the pressure boundary condition 
began 4% aft of the leading edge. Referring to Fig. 22, it can be said 
that although slightly better results were obtained near the leading 
edge for the first ease, that the airfoils designed were fairly insensitive 
to the initial section. 

Tk* Rffeet of the Direct-Inver se Junction 

Proximity to ; l&e R eading Edge 


Since experience with the method has shown that the closer the 
inverse boundary condition is applied to the leading edge, the longer 
it takes for the solution to converge, it was of interest to learn how 
the location of the direct-inverse interface affected the final design an 
the resulting pressure distributions. This study was accomplished with 
the aid of the previously discussed Wing-B case, which began at 2.5 /e 
chord, and an inviscid design of Wing-B also with NACA 0006 sections 
as the initial geometrv. With the second case, the design was begun at 
5% chord from the leading edge, and the input pressures were obtained 
from an inviscid analysis of Wing-B. Since the pressure distributions 
were consistent in both of these cases, the fact that one was a viscous 
design and the other an inviscid design is not important. 

Some representative samples of the resulting section shapes for the 
second case are shown in Fig. 23. The resulting wings were analysed 
‘ under the same conditions that the original input pressure distribu- 
tions were obtained. Representative samples of the resulting pressure 
distributions are compared to their respective target distributions in 
Figs. 24-25. As can be seen, the wing whose design began 2.5% aft 
captured the suction peek at the leading edge, while the other case, 
which began at 5% aft of the leading edge, did not. 

When designing close (less than 5%)' to the leading edge, the 
solution sometimes began to slightly diverge or ceased converging. 
Usually the design could be converged to the point where there was 
only a maximum change in the surface of .l-.2% chord. This was more 
a problem on the fine grid than on the medium. If it was necessary 
to converge it further, the beginning of the design region was moved 
aft. This is an important observation, for if it is necessary to begin the 
design dose to the leading edge to properly determine the shape of the 
nose, a successful design may be accomplished by beginning the design 
as close to the leading edge as desired or is possible, then moving the 
beginning of the design region aft as the solution approaches the last 
stages of convergence. This method not only frees the designer from 
the task of choosing the correct leading edge shape, but it should also 
tVi# ronvetffenee of the desien considerably. 


Fixed Trailine Edge Design 

This case was investigated to verify that a fixed trailing edge 
design could be accomplished with the present version of the code. 
The case chosen utilised Lockheed Wing-A at a Mach number of .8 
and an angle of attack of 2*. A NACA 0012 section was used as the 
initial geometry from 30% to 70% semispan, while the remaining part of 
the wing used the original supercritical sections. The inverse boundary 
condition was enforced from 5% to 80% chord. The airfoil aft of 80% 
chord was fixed so that it maintained the NACA 0012 trailing edge 
shape. The input pressures were obtained through a medium grid 
inviscid analysis of the wing with the original supercritical sections 
used throughout. Furthermore, to provide for a smooth transition at 
the aft direct-inverse junction, the displacements were smoothed in the 
chordwise direction. The type II-2 design method was used in this case. 

The resulting section shapes are revealed in Fig. 26 . The target 
airfoil section would actually be the first 80% of the supercritical section 
and the last 20% of the NACA 0012 section. Surprisingly, even with 
the aft portion of the wing fixed, the designed sections came quite 
close to matching the original Wing-A profiles at the 30% and 50% 
semispan locations. At the 70% semispan location, the designed section 
as compared to the original Wing-A section is much thicker on top and 
thinner on the bottom leading to more cambered profile. This shape 
is probably due to the interaction of the geometric constraints and 
the required design pressures. The shock strength of the input C T 
distribution does become quite large at this location and it appears 
that the section may have had to become more cambered to account 
for this. Or, the increased camber may have been needed to provide the 
necessary lift required by the inverse boundary condition. The pressure 
distributions obained from an inviscid anlaysis of the resulting shapes 


are compared with those produced b^the ongmal wmg ^ ^ design 

the NACA 0012 secuons in Fig. - Th *6 ^ jj aCA 0 012 

pressure distributions are a combination of % & ^ndary 

pressure distributions. It is also interes ng Oial . . thf 

*Hoch near the aft of course, was 

constraints of this P I0 ^^\^ t i "Lible to fix the aft region of the 
'.X idge was used, better result, would 

surely folio*. 

Pressure Bilt ribution CompatibillD: 

, u. .h,t a designer might not readily have 
Since U was though compaL ible with the design 

available an input P" ssu " of intere5t to discover the effect of 

freestream Mach numb * r ' 1 her using a pressure distribution 

designing a wing at one ■ Bt a different Mach number, 

obtained from an analyisis h 8 hout this portion of the study. 

:^tery= t oot * *, **- -*■ ~ 

begun at 107c aft of the leading edg involved a fine 

Two separate number of .2 using a pressure 

design at a nearly incompr target *u a Mach number 

» r*' J ,rv:: 

of .1. As can be seen from Fig. - , . -.a the j.p Prandtl- 

at the higher Mach number. This rs rn agreement with 

Glauert similarity rule* 

(35) 



which states that the C r will 

thickness, r is reduce as te ^a^ ^ th&t ,. 54 % decrease in 

flow. For this case, Eq. ( ) f prc ssure distribution at 

thickness would be ° codc f ot this 3-D case produced 

the higher Mach ™ mbtt , T g % lhinner th „ the NACA 0012 

a section which was on the average i.vs 

section. . . * grid design at a Mach number 

The second case involved a medium grid e g of >80 . 

of .85 using a pressure distribution o ^ ^ again thinner than 

Referring to Fig- 29, the section s ap * requ ired a sudden thinning 

the initial section. The top surface. though, analysing this 

of the surface at the shock location. Surpris g ; P with 

the target everywhere excep boundary 

“ f - - 111 

pressure distribution to eliminate the necessity of these dips. 

Grid Refinemen t Effects 

■ Si„« a. «»• £ -fiK ZZX 

grid is about an eighth of that teB l pressures. In 

to try to design on the medmm S^-g fi^gnd P & medi 

order to assess the practicality ■ Th ease waJ performed 

grid using fine grid pressure, “ atli of two degrees. The 

using a Mach number of .8, and on g uJtd u the initial, as 

original supercritical sections ^or^V/^g ^ p . g 31 xhe only' 

well as, the target sections. target was near the middle of 

place where the designs came of the designed 

the wing. A slight wave j^ to the smearing of the shock 

vJzfiST* «*? ”° 5t ptow>,y nectssiuie<1 the 

decambering of the sections at the wing tip. 

No attempt was made to match the C.'s of the fine grid and 
medium grid analyses * or angle of aU^ 

S Zl~V to Mtef £ ZL of g the wing to closely match 

of the design stations. It also shows that increasing the angle of attack 


, • ,, hBVt produced closer matching C,’s and hence perhaps 

irfouni. to' use the fine grid to properly resolve the 
correct airfoil sections. 

Conclusions 

hr a. — ^5.*tol Ml M.. 

— n I— - 

be performed on a fine grid. 
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Fig. 3 The effect of relofting on the design 
in the initial stages of convergence 



Summation of the Residual Terms 




Fig. 5 The time history of the residual and the terms composing it 
at the uppper trailing edge for a typical divergent solution. (The 
design region extends from 30-70% semispan) 
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Figs. 6a-6e Comparison of airfoil sections designed using the four 
spanwise oscillation remedies. 
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Fig. 6d 




Fig. 6b 




Fig. 6c 


Fig. 7 The effect of the four spanwise oscillation remedies on the 
average percent error between the target and the design sections 
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Summation of the Residual Terms 





Fig. 1H 



Fig. lib 




Fig. lid 


Figs, 1 la-1 Id Comparison of a pressure distribution obtained from 
an analysis of a skewed grid with one obtained from an analysis o 
an unskewed grid 




Fig. 12c 

Figs. 12a-12c Comparison of the sectional shapes designed using 
different viscous interaction assumptions with a pressure distibu- 
tion obtained from a fully viscous analysis of Lockheed Wing-A at 
a Re = 24x10* , M = .8, and an a = 2 ° 
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Fig. 13a 

Figs. 13a- 13c Comparison of pressure distributions obtained from 
a viscous analysis of the Lockheed Wing-A using three drfferent 
viscous interaction assumptions with a Mach - .8, a , 

Re = 24x10* 
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Fig. 15a 


Fig. 13c 
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Fig. 16 Grid generated about VVing-C with an incompatible root 
section and fuselage cross section 



Fig. 17 Surface velocity vectors for Lockheed Wing-A on the upper 
and lower surface 
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Fig. 19b 
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Fig. 19c 



Fig. 18 Surface velocity vectors for Lockheed Wing-C on the upper 
and lower surface 


Figs. 19a- 19c Comparison of the designed sections with the targets 
and the initial sections for a fine grid case using Lockheed Wing-A 



Figs. 20a-20c Comparison of the designed sections with the targets 
and the initial sections for a fine grid case using Lockheed Wing-B 
and a design region beginning at 2.5% aft of the leading edge. 


H 





0-VO 









.g 000 
o 
o 

V) -0.07 


0.0 6, « o.3 0.4 oT o.e 03 o.. m vo 
Chord Fraction ix/c) 





/ 

Figs. 22a-22c Comparison of sections designed using two difforent 
initial sections 
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Fig. 21b 











Fig. 23b 



Figs. 23a-23c Comparison of the designed sections with the targets 
and the initial sections for a fine grid case using Lockheed Wing-B 
and a design region beginning at 5.0% aft of the leading edge. 



Figs. 24a-24c Comparison of the pressure distributions obtained 
from an analysis of the Wing-B design, which had a design region 
which began 2.5% aft of the leading edge, with the target pressure 
distributions 



Fig. 24b 





Figs. 25a-25c Comparison of the pressure distributions obtained 
from an analysis of the Wing-B design, which had a design region 
which began 5.0% aft of the leading edge, with the target pressure 
distributions 
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Figs 26a-26c Comparison of the sectional shapes designed with 
a 1 fixed trailing edge region with a NACA0012 and the original 
Lockheed Wing-A section 



Figs. 27a-27d Comparison of the pressure a fixed 

an inviscid analysis of the . sect '°^‘ ^^CAOOU and the original 
trailing edge region with those for a N ACAO . 

Lockheed Wing-A section. (M -8> Q - > 
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Fig. 29c 



Figs. 29a-29c Comparison of the sections designed at a Mach 
number of .85, using input pressure distributions obtained from 
an analysis at a Mach number of .8, with the original supercritical 
sections 



Fig. 30a 

30a-30c Comparison of the pressure distributions obtained 
, an analysis at a M = .85 of the Lockheed W.ng-A which iwm 
ened at a M = -85 using input pressure distributions obtained 
f an analysis at a Mach number of .8. with the input pressure 
ributions 


Figs. 31a~31d Comparison of the sections obtained from a medium 
grid design at at M = .8 using input pressure distributions obtained 
from a fine grid analysis of Lockheed Wing-A with the target 
sections 
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INVERSE WING DESIGN IN TRANSONIC FLOW INCLUDING VISCOUS INTERACTION* 
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|| SUMMARY 

|j li- 

jgp; Several inverse methods have been compared and initial results indicate 

a|f that differences in results are primarily due to coordinate systems and 
fuselage representations and not to design procedures. Further, results 
from a direct- inverse method that includes three-dimensional wing boundary- 
layer effects, wake curvature, and wake displacement are presented. These 
“ results show that boundary-layer displacements must be included in the 
pT design process for accurate results. 


1 


^ INTRODUCTION 

Over the past several years, a variety of transonic wing design methods 
and computer codes (refs. 1-5) have been developed. In general, these 
methods solve the full potential flow equation and utilize the inverse 
approach in that pressure distributions are specified over all or part of 
the wing surface. Several include some of the effects of viscous 
interaction via strip boundary-layer calculations (ref. 1) or two- 
dimensional computations that include a correction for three-dimensional 
viscous effects (ref. 3). However, none of these methods includes a true 
three-dimensional boundary— layer calculation or the effects due to wake 
curvature, etc., which might have important effects on computed wing 
designs . In addition, they differ in the number and spacing of grid points, 
the design approach, the treatment of fuselage effects, and the control of 
ppP trailing— edge thickness . Obviously whether or not these formulation 
differences significantly affect design results is of interest. 

swl---- Currently, the design version of TAWFIVE (refs. 6-7), termed TAW5D (ref. 
4), is being extended to include three-dimensional boundary -layer and wake 
viscous interaction effects and is being used to study various leading— edge 
relofting/trailing— edge control design procedures. As part of this study, 
it was believed that it would be interesting to investigate the consequences 
of differences in both numerical and physical formulations on the design 
process and resultant wing designs. Thus, this paper will present initial 
results of two ongoing studies. The first part will compare several inverse 

*This work was supported by NASA Grant NSG 1-619. 









ORIGINAL FACiil ’3 

OF POOR QUALITY 


design methods and their results, while the second portion will discuss the 
influence of viscous interaction on transonic wing design. 

INVERSE METHOD COMPARISON STUDIES 

The RAE Wing Body 'A' configuration (ref. 8) at a freestream Mach 
number of 0.8 and angle of attack of 2 degrees was selected as the test case 
for the comparison studies. The wing for this configuration has an aspect 
ratio of 5.5, a leading-edge sweep of 36.7 degrees, and a taper ratio of 
0.375, is untwisted, and is composed of RAE 101 symmetrical airfoil 
sections. Three different inverse design methods were selected for the 
comparison, the direct- inverse curvilinear coordinate system TAW5D code 
(ref. 4), the stretched Cartesian grid direct- Inverse ZEBRA method (ref. 
2-3), and the inverse predictor-corrector FL030DC approach (ref. 5); and 
their characteristics and features are listed on Table I. 

In order to avoid the complexities associated with various viscous 
interaction schemes, it was decided to limit this comparison study to 
inviscid flow; and, since it was believed that one of the primary usages of 
design codes would be to modify only portions of wings, it was decided to 
design only between 30 and 70 percent span. The target pressure 
distribution for the design zone was obtained from an inviscid analysis by 
the TAW5D code (essentially TAWFIVE, ref. 7), which indicated that the 
flowfield at the selected conditions was slightly supercritical and that the 
wing lift coefficient was 0.210. In addition, the starting airfoil shapes 
were the correct 9% thick sections from root to 30% span, linearly thinning 
down to a 6% thick symmetrical section at 50% span and back to 9% at 70% 
span, followed by the correct sections on the outboard portions of the wing. 

For the design studies, TAW5D was operated in the span lofting mode in 
which pressures were only specified at 30, 50, and 70% span. Under this 
procedure, airfoils were only Inversely designed at these stations; and 
after each design update, in between sections were obtained by linear 
spanwise lofting. In all cases, the flow at these in between stations was 
computed in the direct -analysis mode. On the other hand, in the ZEBRA 
method, pressures were specified at each spanwise station from 30 thru 70%; 
and In the predictor-corrector, FL030DC the pressure was specified and an 
airfoil section designed only at the 50% span location, with linear span 
lofting to 30% and 70% respectively. In all cases, leading-edge relofting 
options were selected in order to force the designs to have the proper 
trailing-edge thicknesses . 


PROBLEMS 

In setting up the test cases, several interesting problems were 
encountered. First, analysis computations of the RAE 'A' wing/body 
configuration by the ZEBRA and TAW5D codes yielded slightly different 
pressure distributions; and, in order to minimize these differences, the 
angle of attack used in ZEBRA was decreased to 1.8 degrees so as to match 


I wine CL predicted by TAW5D . The corresponding pressure distributions 
sho^n on figure 1; and since both methods solve the same equation, the 
-ust be due to differences in grid, fuselage, and boundary 
Edition treatments. Near the root, ZEBRA predicts a greater . 

that the flow is more accelerated on the upper surface, while 
Aboard the leading-edge grid clustering inherent in TAW5D results in 
Ktter resolution of the leading-edge region and minimum pressure P** k - 
t l the trailing edge, where the ZEBRA coordinate system is actually finer , 
Sere are a?so some variations in the predicted pressures. However, between 
Jfo'and 70% span the two methods are in reasonable agreement, and meaning 
Resign studies for this region should be possible. 

The second problem was that FL030DC could only handle for this case an 
finite cylinder fuselage; and, thus, TAW5D and ZEBRA were modified to 
^ve as an option an infinite fuselage as well as a finite one. F gu 
L' oa res at the 50% span station on the RAE configuration the pressu 
KisLibutions calculated by TAW5D associated with these two fuselages, and 
KtZl be seen that the effect is only a slight shift in the V ressnre 
Coefficient level. This trend was true at all span stations and overall 
# intr and section lift coefficients were essentially identical. 

fgevLSeless es a result of these differences two sets of target pressures 
*?fr>r the design region were generated, one for the finite wing/bo y 
l^onf iguration and'one for the infinite cylinder fuselage; and these were 
|used as input into the appropriate versions of the codes. 

r 

Results and comparisons 

Hr" Figures 3-5 show results obtained at the design stations using the 
ItaV 5D method In this case, each section was designed from 10% chord to th 
Ittailing edge and leading-edge relofting was utilized to force trailmg-edge 
TVlosure 6 However the actual ordinate of the trailing edge was not 
fed. "”; be seen, tha starring profiles were a linear variation 

ffrom the correct section at 30% and 70% span down to a thin 

^section at mid-span. While the 30 and 70% stations started with the correct 
thanes they were design stations and could and did change during the 
g^UtSn 7 However, S .s shown on the figures, all three sections converged 
|“to the target shapes; and results for the finite fuselage and infinite 
cases were indistinguishable . 

PS*- Results were also obtained with the ZEBRA code for both the infinite 
land finite body cases and by the FL030DC code for the infinite ^nder 
Tfuselage using the appropriate pressure inputs. Figures &- 8 compare the 
^designed sectional shapes obtained by the three codes for the infi ^® 
Ifuselage It should be noted that the ZEBRA results were well converged 
if avi^g maximum ordinate changes of less than IE- 6 of chord when computations 
Were terminated. Also, it can be seen that the FL030DC and TAW5D resul 
►(denoted as CAMPBELL and TAWFIVE on the figures) are virtually identical, 
hven though the methods used entirely different design procedures. 
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At the 30% span station, the lower surface profile predicted by ZEBRA 
is in agreement with the other methods, but on the upper surface it is 
considerably different. Examination of the pressure profiles on figure 1 
indicate that at 30% TAW5D and ZEBRA analysis results agree on the lower 
surface but disagree on the upper. Consequently, when the TAW5D pressures 
are used as design input to ZEBRA, it is not surprising that a slightly 
different airfoil section resulted. At 50%, figure 7, where analysis 
results are in better agreement, particularly on the upper surface, the 
three methods predict virtually identical upper surfaces although the ZEBRA 
lower surface profile is slightly thicker; and at 70% span the ZEBRA 
prediction is again slightly thicker. (Similar differences between TAW5D 
and ZEBRA were obtained for the finite fuselage case.) Since TAV5D and 
ZEBRA use similar design procedures and TAW5D and FL030DC have similar grids 
and body representations , it can be concluded that the differences in 
profile shapes portrayed in figures 6-8 are primarily due to coordinate 
system and fuselage representations. 

In order to see infinite versus finite fuselage effects, the infinite 
cylinder fuselage wing pressures were used as input into both the infinite 
cylinder and wing/body versions of TAW5D; and a typical result is shown on 
figure 9. Here the infinite cylinder result is the "correct 11 profile; and 
as can be seen, the finite fuselage result is thinner and significantly 
different near the trailing edge. In fact, at the 30 and 70% stations, the 
upper and lower surfaces criss-crossed before coming together to satisfy 
trailing-edge closure. It is believed that this result demonstrates an 
important effect often encountered in inverse design, i.e., when a pressure 
distribution that is somehow incompatible with either physical reality or 
the computational model (in this case the fuselage representation) is used 
as input, the effect is almost always observed as either unrealistic 
profiles near the trailing edge or in the inability of the design process to 
satisfy the design input pressures near the trailing edge or both. In many 
cases, the "problem" can be solved by slight adjustments in the specified 
pressure distribution. 

Now even though figures 6-8 show that the methods predicted different 
profiles, the significance of these differences can only be determined by an 
analysis of the designed wings and a comparison of the analysis results with 
the desired targets. Since TAW5D had previously been shown to be self 
consistent (ref. 4) and since the wing designed by TAW5D, fig. 3-5, had the 
correct airfoil sections, no analysis results for the TAV5D design are 
presented. However, figures 10-14 compare the target pressure distributions 
with analysis results by both TAW5D and ZEBRA for the wing designed by 
ZEBRA, which had different profile sections in the design region. First, it 
should be noted that in the design region, figures 11-13, the ZEBRA analysis 
agrees with the target pressure values for the inverse design zone, which 
extends from 0.1 chord to the trailing edge. This agreement indicates that 
the ZEBRA method did indeed satisfy the desired pressure boundary 
conditions. Second, due to inherent grid clustering near the leading edge, 
the TAW5D analysis of the ZEBRA design probably gives better resolution in 
the leading-edge region; and, finally, if it is assumed that the TAW5D 
analysis is the "most accurate" of the methods due to its fuselage and 
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are shown on figures 15-17, and as can be seen the initial sections linearly 
varied from the correct aft-cambered profile at 20% span to a conventional 
non-aft cambered section at mid-span back to the correct aft-cambered 
section at 80% span. Here, the inverse design procedure started at 0.1 
chord; and the initial leading edge at each design station was thinner than 
the target shape. As shown on the figures, the target sections and designed 
sections are in excellent agreement, particularly considering the extensive 
curve fits and interpolations involved in the design and viscous interaction 
procedures . 

For the second test, the initial sections consisted of the correct 
profiles inboard from the root to 20% and outboard from 80% to the wing tip. 
However, as shown on figures 18-20, from 30% span through 70% span the initial 
sections were NACA 0012 airfoils; and linear lofting was used between 20 and 
30% and 70 and 80%. In this case the inverse design procedure started at 
0.04 chord, and the initial leading edge at each design station was thicker 
than the target section. As can be seen, the final designed sections are in 
excellent agreement with the target shapes, particularly in the leading-edge 
and cove regions. 

It should be noted that in both of these cases, the section and wing 
lift coefficients and the section pressure distributions were essentially 
identical to the target values. Based upon these results, it is believed 
that the present viscous inverse design procedure can yield correct target 

iles independent of initial airfoil section shapes . 

BOUNDARY- LAYER AND WAKE EFFECTS 

Studies conducted under the present program have indicated that design 
including full viscous interaction effects is more computationally intensive 
and that convergence is slower. Consequently, it was decided to compare the 
full viscous interaction design results with those obtained including 
viscous boundary-layer interaction but excluding wake effects and with those 
obtained assuming inviscid flow. For each case, the input pressure 
distributions were identical and corresponded to those predicted by a full 
viscous analysis of the Lockheed Wing A wing/body since those should be the 
closest to reality. The starting section profiles were those shown on 
figures 15-17, and the design region was from 30 to 70% span. As before, 
span relofting and leading-edge relofting were both used in all three cases. 

The final section profiles resulting from these computations are shown 
on figures 21-23, and at all design stations the sections obtained by 
ignoring wake effects are very close but slightly thicker than those 
corresponding to the full viscous case. Further, while the inviscid case 
profile is very close to the others at 50% span, they are significantly 
different from those including viscous effects at 30 and 70% span. The 
results at 50% are not surprising since at that station the boundary layer 
is relatively thin over much of the surface and the design is strongly 
influenced by the viscous pressure boundary conditions at 30 and 70% span. 
However, the cove region is not well predicted; and, as can be seen on 
figure 22, the upper surface inviscid profile here is thinner than the full 
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CONCLUDING REMARKS 

rec.lJc fT/T Saveral inverse methods have been compared and initial 

systems and fusel a differenCe f in results are primarily due to coordinate 
results from an i ^ re P res <-ntations and not to design procedures. Also, 

lavei effects waW that includes three dimensional wing boundary- 

layer effects, wake curvature, and wake displacement have been presented 

These results show that boundary-layer displacements must be included in the 
design process for accurate results. 
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TABLE I. -- CHARACTERISTICS OF INVERSE METHODS 


Method 

TAW5D 

ZEBRA 

FL030DC 

Coordinate System 

Body Fitted 

Stretched Cartesian 

Body Fitted 

Boundary Conditions 

On Surface 

At Z - 0 

On Surface 

Fuselage 

General Shape 

Axisymmetric Body 
Approx, by Source/Sinks 

Infinite 

Cylinder 

Design Method 

! 

Direct-Inverse 

Direct- Inverse 

Predictor- 

Corrector 

Grid 

160x24x32 

90x30x30 

160x24x32 

Points on Airfoil 
Section 

105 with LE 
Clustering 

100 almost equally 
spaced 

105 with LE 
Clustering 

Number of Span 
Stations 

21 

21 

21 


TABLE II. -- RESULTS OF ANALYSIS OF DESIGNED WINGS 



Target 

Full Viscous Design 

No Wake Design 

Inviscid Design 

cl at 50% 

.514 

.509 

.506 

.427 

Wing CL 

.483 

.478 

.477 

.419 
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Figure 3. Comparison of section designed by TAW5D at 30 percent span for 
RAE wing body 'A' with initial and target sections. 




Figure 4. Comparison of section designed by TAW5D at 50 percent span for 
RAE wing body 'A' with initial and target sections. 
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5. Comparison of section designed by TAW5D at 70 percent span for 
RAE wing body 'A' with initial and target sections. 
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Figure 6. Comparison of sections designed by different methods at 30 
percent span for RAE wing 'A' with infinite fuselage. 
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re 7. Comparison of sections designed by different methods at 50 
percent span for RAE wing 'A' with infinite fuselage. 



Figure 8. Comparison of sections designed by different methods at 70 
percent span for RAE wing 'A' with infinite fuselage. 
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Comparison of sections designed by finite and infinite fuselage 
versions of TAW5D using infinite fuselage wing pressures as input 
in both cases (RAE wing 'A'). 
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Comparison at 10 percent span of target values with pressures 
obtained by analyses of the wing designed by ZEBRA 
(RAE wing body 'A', Mach — 0.8, AOA — 2 degrees). 
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Figure 11. Comparison at 30 percent span of target values with pressures 
6 obtained by analyses of the wing designed by ZEBRA 

(RAE wing body 'A', Mach - 0.8, AOA - 2 degrees). 
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Comparison at 50 parent span of target "“ h P« ssures 

obtained by analyses of the ving designed by ZEBRA 
(RAE ving body 'A' , Math - 0 . 8 , AOA - 2 degrees). 
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13. Comparison at 70 percent span of target values with pressures 
obtained by analyses of the wing designed by ZEBRA 
(RAE wing body 'A' , Mach - 0.8, AOA - 2 degrees) . 


- 0.4 4 / 


90 PERCENT SPAN 


CP DISTRIBUTIONS 
TARGET 
TAWRVE 
"ZEBRA"” 


TARGET Cy = 0.222 
TAWFIVE C l « 0.217 
ZEBRA C, = 0.212 


0.3 0.4 0.6 0.8 0.7 

CHORD FRACTION, X/C 


Figure 14. Comparison at 90 percent span of target values with pressures 
obtained by analyses of the wing designed by ZEBRA 
(RAE wing body 'A' , Mach - 0.8, AOA - 2 degrees). 
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15. Comparison of section designed by TAW5D at 30 percent span for 
Lockheed Wing A wing body with target and first type of initial 
section. 
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Figure 16. Comparison of section designed by TAW5D at 50 percent span for 
Lockheed Wing A wing body with target and first type of initial 
section. 
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Figure 17. Comparison of section designed by TAU5D at 70 percent span for 
Lockheed Wing A wing body with target and first type of initial 
section. 
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Figure 18. Comparison of section designed by TAW5D at 30 percent span for 

Lockheed Wing A wing body with target and second type of initial 
section. 
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19. Comparison of section designed by TAW5D at 50 percent span for 

Lockheed Wing A wing body with target and second type of initial 
section. 
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20. Comparison of section designed by TAW5D at 70 percent span for 

Lockheed Wing A wing body with target and second type of initial 
section. 
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are 21. Comparison of sections designed at 30 percent span using 
different viscous interaction assumptions for Lockheed 
Wing 'A' wing body. 
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Figure 22. Comparison of sections designed at 50 percent span using 
different viscous interaction assumptions for Lockheed 
Wing 'A' wing body. 
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Comparison of sections designed at 70 percent span using 
different viscous interaction assumptions for Lockheed 
Wing 'A' wing body. 
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Comparison of pressures at 10 percent span obtained by viscous 
analyses of the wings designed using dl ff® re ^ Z_ os 
interaction assumptions (Lockheed Wing , a ’ 

A0A - 2 degrees, Reynolds No. - 24 million). 
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25. Comparison of pressures at 30 percent span obtained by viscous 
analyses of the wings designed using different viscous 
interaction assumptions (Lockheed Wing 'A', Mach - 0.8, 

AOA - 2 degrees, Reynolds No. - 24 million). 



Figure 26. Comparison of pressures at 50 percent span obtained by viscous 
analyses of the wings designed using different viscous 
interaction assumptions (Lockheed Wing 'A', Mach - 0.8, 

AOA - 2 degrees, Reynolds No. - 24 million). 
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Figure 27. Comparison of pressures at 70 percent span obtained by viscous 
analyses of the wings designed using different viscous 
interaction assumptions (Lockheed Wing 'A', Mach - 0.8, 

AOA - 2 degrees, Reynolds No. - 24 million). 
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Figure 28 . Comparison of pressures at 90 percent span obtained by viscous 
analyses of the wings designed using different viscous 
interaction assumptions (Lockheed Wing # A # , Mach - 0.8, 

AOA - 2 degrees, Reynolds No. - 24 million). 
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abstract 

Inverse Transonic Wing Design Using Finite -Volume 
Methods in Curvilinear Coordinates (May 1987) 

Thomas Anthony Gaily, B.S., Texas A&M University 
Chairman of Advisory Committee: Dr. Leland A. Carlson 

An inverse wing design method has been developed around an existing 
transonic wing analysis code. The original analysis code, TAWFIVE, has 
as its core the numerical potential flow solver, FL030, developed by 
Jameson and Caughey. Features of the analysis code include a finite - 
volume formulation; wing and fuselage fitted, curvilinear grid mesh; and 
a viscous boundary layer correction that also accounts for viscous wake 
thickness and curvature. The development of the inverse methods as an 
extension of previous methods existing for design in Cartesian coordi- 
nates is presented. Demonstrative results are shown for inviscid wing 
design in both sub-critical and super-critical flow regimes. The test 
cases selected also demonstrate the versatility of the design method in 
the designing of an entire wing or any discontinuous sections of a wing. 
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INTRODUCTION 


In recent years the importance of transonic flight to both military 
and commercial aircraft and the development of specialized transonic 
wings for several flight research experiments have prompted significant 
efforts to develop accurate and reliable computational methods for the 
analysis and design of transonic wings. Hany methods of solution have 
been developed, but those which have shown the most promise due to their 
computational efficiency and engineering accuracy have been based upon 
the full potential flow equations in either their conservative or non- 
conservative form*’ 3 . When these potential flow codes have been coupled 
with accurate viscous boundary layer routines, they have had great suc- 
cess in predicting such complex flow phenomena as wing-body interac- 
tions, three-dimensional shock wave formations, and weak viscous inter- 

actions . 

The TAWFIVE 4 FORTRAN code in particular has proven to be an excel- 
lent and reliable analysis tool. This analysis code is based upon the 
FLO 30 finite volume potential flow method that was developed by Jameson 
and Caughey 3 . Among the features of FL030 are its fully conservative 
formulation and its three-dimensional curvilinear grid. The latter can 
be fit around any general combination of fuselage shape and wing plan- 
form. In addition, TAWFIVE gives the user the option of including the 
viscous effects associated with both a wing surface boundary layer and 
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viscous wake having both finite thickness and curvature. 

The purpose of the research outlined in this thesis has been to 
develop an inverse wing design method that is based on the existing TAW- 
FIVE analysis code and is compatible with the existing computational 
methods and program structure of that code. The particular inverse 
method used extends previously developed design methods^' ^ developed for 
orthogonal grids to the more generalized curvilinear grid system of TAW- 
FIVE, while also providing for greater design flexibility and versatil- 
ity for engineering applications. These last goals were achieved by the 
inclusion of user options for designing either the entire wing or only 
discontinuous wing segments as shown in Figure 1. The availability of 
this option is useful to engineers who are typically faced with desig- 
ning around regions where the wing geometry may be fixed by constraints 
other than aerodynamic considerations . 



(a) Part of Upper Surface, 
Lower Surface, or Both. 





Figure 


1. Possible Wing Design 


Situations 
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BACKGROUND 

There are a number of approaches to the solution of viscous tran- 
sonic flows over wing-body combinations. The most exact method, if a 
suitable numerical algorithm is available, would be to solve the com- 
plete set of continuity, momentum, and energy equations with viscous 
effects either included in the governing equations or with a separate 
boundary layer solver; i.e. the Navier- Stokes or Euler equations with 
boundary layer. At present, however, solution methods of this type are 
still under development and are far too demanding of computational speed 
and memory requirements to be used for practical engineering applica- 
tions. However, with the great advances being made in both numerical 
methods and computer capacity, solutions of this type may become common 
in the future . 

As previously mentioned, the solution methods which are primarily 
used today combine the compressible potential flow equations with a vis- 
cous boundary layer solver. The compressible potential flow equations 
are : 

x x y y (1) 

P - [1 + - *X 2 • V - *z 2 >! 7 1 

The derivation of these equations from the continuity and inviscid 
momentum and energy equations can be found in almost any textbook deal- 
ing with compressible subsonic or transonic flow. The use of the invis- 
cid forms of the governing equations as the basis for the potential 
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equations is necessary due to the fact that the potential function is 
not defined for fluid flow in the presence of vorticity. the natural 
by-product of viscous effects. While the potential equations are suit- 
able for most flow regions where viscous effects are negligible, there 
are two situations where viscous effects cannot be ignored. The first 
is the boundary layer over a wing, and the second is across shock waves. 

The wing viscous interaction problen is usually solved by separating 
the flow field into two regions: the inviscid region away from the body 
where the potential flow equations apply and the viscous region on the 
wing surface where a separate set of boundary layer equations can be 
solved. The boundary for the potential solution region is usually off- 
set from the physical surface by a distance called the boundary layer 
displacement thickness . This thickness is calculated within 
ary layer routines based upon the momentum loss and mass flow decrease 
of Che flow within the viscous layer. Since the potential flow 
boundary layer equations are usually solved separately but are dependent 
upon each other, it is necessary to iterate between the two methods in 

order to obtain a final converged solution. 

The second viscous region of concern, shock waves, is accounted for 
by the approximate numerical solution methods used in transonic poten- 
tial flows. An examination of the nonlinear partial differential poten- 
tial flow equations reveals that the nature of the equations change 
fro. elliptic to hyperbolic as the flow accelerates from subsonic to 
supersonic speeds. In numerically solving these equations of mixed 
nature , it is necessary to use a numerical approximation which exhibits 
the proper upstream dependence in the supersonic flow regions. 




In the 
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original paper by Murman and Cole^ suggesting this flow dependent method 
for solving transonic flows, the authors noted that their numerical 
solution technique for supersonic flow regions introduced a numerical 
error term to the inviscid equations which had the effect of an artifi- 
cial viscosity. The result was the appearance of numerical shock waves 
in their solutions which were very similar in nature and location to 
what would be expected of physical shock waves. While this result is 
somewhat serendipitous, it is of great consequence and has become the 
foundation of transonic potential methods. Methods newer than those of 
Murman and Cole do not make use of the same numerical approach but 
always try to introduce a supersonic artificial viscosity of roughly the 
same character. 

In the area of numerical airfoil/wing design using transonic poten- 
tial flow solutions, recent efforts can be categorized as being either 
optimization or inverse methods. The optimization problem can be 
expressed mathematically as finding the local extrema (maxima or minima) 
of a function of many variables subject to a given set of constraints. 

A practical example would be to find the set of airfoil ordinates which 
would produce the minimum transonic wave drag but which is not undesir- 
able thin. 

For a function described by a set of analytic expressions, this 
optimization problem is straight forward and has distinct solutions. 

When applied to airfoil/wing design, the problem is complicated by the 
nonlinear, often discontinuous, and unknown form of the dependent func- 
tion and also by the numerical accuracy and nature of the flow field 
solver to which the optimization routine is coupled. With a proper 
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selection of independent variables and constraints, however, the problem 
becomes manageable and good results can be obtained. 

An excellent example of the optimization approach is the work of 
Cosentino and Holst*. They chose to use a finite number of wing surface 
ordinates as their independent functions and exercised their program by 
minimizing the ratio of drag to lift for various wing planforms and 
flight conditions. However, most airfoils require a large number of 
surface ordinates, on the order of 30 or more, to smoothly define the 
airfoil surface. For a typical wing which may be defined by 10 or more 
spanwise airfoil ' sections , this approach could possibly involve 300 
independent variables. As mentioned earlier, the relationship between a 
change in a surface ordinate and the resulting change in the property 
being minimized is not known a priori but must be calculated from the 
flow solution itself. The relationship can be determined by displacing 
each of the surface ordinates individually and evaluating the response 
of the dependent function. If all the relationships were linear in 
nature, a minimum of N iterations would be necessary to minimize a func- 
tion of N variables; and more than 300 converged flow fi'eld solutions 
would have to be obtained in order to minimize the proscribed optimiza- 
tion function. Obviously, such a procedure would require an unaccept- 
able amount of computational time. In practice, however, the optimiza- 
tion routine attempts to drive the solution to the minimum before the 
system has been completely defined; and it is quite probable that the 
minimum can be reached before N iterations. Nevertheless, a large num- 
ber of iterations can still be expected for unconstrained problems. 

To reduce the size of the task, Cosentino and Holst only varied a 
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few surface ordinates over a restricted region of the upper surface and 
depended upon spline fits through these points to provide a smooth sur- 
face. Of course, by minimizing the number of surface ordinates, the 
problem is constrained to only having solutions within the family of 
spline curves which can be fit through the points. 

Another recent work by Davis^ reduces the airfoil surface to a 
series of parametric curves that have predetermined functional relation- 
ships to certain flow phenomena. It only remains for the optimization 
routine to determine the functional relationship between the specific 
flow phenomena and the desired optimization function. For example, if a 
family of curves for a section of say the upper surface has been found 
to have a particular effect on the location of shock waves, then if a 
location for a shock wave can be found which maximizes the lift to drag 
ratio for the airfoil, the surface which corresponds to this shock loca- 
tion is proscribed by the function relating the two. In essence this 
approach is another means of constraining the problem to manageable lev- 
els as was done by Cosentino and Holst, but it is accomplished in a more 
elaborate manner and requires and makes use of known surface to flow 
relationships . 

The other set of design techniques are called inverse methods. 

Their distinguishing feature is that they require the specification of 
desired velocity or pressure distributions over an airfoil or wing and 
calculate the corresponding surface shapes. The usefulness of inverse 
methods lies in the fact that experienced airfoil designers will usually 
have an idea of what velocity /pressure distributions will produce desir- 
able design properties. An example of such design would be modifying 

i* 
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the pressure distribution on an existing wing section so that an unde- 
sirable shock could be moved or even eliminated. Such design cases will 
typically occur when an airfoil has been optimized for a given flight 
condition but is also expected to perform adequately in off design con- 
ditions . 

Within the range of inverse methods, there exist two distinct sub- 
groups. For one group of methods, the transonic flow field solution in 
used as a "black box" that calculates the velocity/pressure distribution 
for a given geometry. The inverse procedure then generates a modified 
geometry based upon the difference between the calculated and desired 
flow profiles. The NYU code of Bauer, Garabedian. and McFadden 10 and 
the more recent method of Takanashi 11 are two such inverse techniques. 

In a sense , these types of inverse methods are related to the optimiza 
tion methods with the minimized function being the difference between 
the actual and desired velocity/pressure profile. These inverse methods 
differ from the optimization techniques because, instead of allowing the 
functional relationship between shape profile and flow field to be 
determined numerically from the convergence history, the functional 
relationship has been determined beforehand based upon the governing 
flow equations . 

The other method of inverse wing design, called the direct- inverse 
method, uses the desired pressure distribution for an airfoil or wing 
surface as a boundary condition which must be incorporated into the flow 
field solution. Thus, rather than satisfying the usual flow tangency 
condition on the body surface, the flow instead must satisfy the desired 
pressure distribution in the inverse regions. 




From the flow direction 
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on the pressure boundary or from the mass flux through it, the location 
of a physical surface satisfying the desired pressure profile can be 
calculated. Since the equations will generally only be exact when the 
pressure boundary and physical boundary coincide, these methods also 
involve some iteration on the surface geometry. 

As mentioned, the primary task of this research effort has been to 
modify the existing finite volume transonic potential flow analysis code 
known as TAWFIVE for use in wing design. In selecting a design method 
to apply to this existing code, it was decided use the direct- inverse 
method as developed by Carlson^ * ® ^ ^ . Carlson's previous inverse 
design work has been with both airfoil design^* ^ and, more recently, 
with wing design ^ , The reasons for selecting this method include the 
relative simplicity of the method along with its record for successful 
use with a variety of different numerical algorithms. In addition, the 
previous applications of this method have been in Cartesian or sheared 
Cartesian coordinate systems. Thus, by adapting this method to the cur- 
vilinear coordinate system of TAWFIVE, original work and results would 
be obtained. 


H o 


11 


WING ANALYSIS METHODS 


Potential Flow Solver 

The inviscid potential analysis of TAWFIVE is performed by the pro- 
gram FL030 developed by Caughey and Jameson 3 - 14 . For a comple-e 
description of the FL030 code and its theoretical basis the reader is 
referred to Caughey and Jameson's papers and some earlier developmental 
work by Jameson 15 ' 16 . A brief description is presented here to provide 
for completeness and to provide a background for the inverse design 

developments which will be discussed in detail. 

PL030 solves Che full potential equation in conservative form which 
when transformed from Cartesian coordinates, Eq. (1). to generalized 

curvilinear coordinates is : 

(phU)^ + (phV ) ^ + (phW)j- - 0 

where p is the local density; U, V, and W are the components of the 
contravariant velocity vector; h is the Jacobian of the transformation 
from Cartesian coordinates; and subscripts denote differentiation with 
respect to the curvilinear coordinates r \ , and f . The contravariant 
velocities are related to the physical velocities and the derivatives of 

the potential function by: 



( 2 ) 


and H is the transformation matrix defined by. 
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The local density can be obtained from isentropic relations as: 

1 

p - [1 + * u 2 - v 2 - v 2 )] 7 ( 4 ) 

The numerical approach used in FL030 is a finite volume technique. 

To understand this approach , consider the simple two dimensional case 
represented by the grid system shown in Figure 2. The dashed cube shown 
in the figure indicates the area element under consideration. The flux 
of fluid through side a-b can be approximated by the average of the 
fluxes at point a and b with similar results for the side c-d. The net 
flux in the x direction for the elemental area centered at grid point 
i t j is then: 

(phU)£ - [<phU a + phU b ) - (phU c + phU d ) ] / 2Af 
or in the notation of Caughey and Jameson, 

(phU) ? - /i n 5 f (pU) 

where p indicates averaging and S indicates differentiation in the 
indicated directions which are defined as follows (allowing A£ -A^-Af-l) : 

( 5 £ u >i,j,k “ ( u i+h , j ,k * U i-H,j,k) 

t k " < u i+H,j ,k + ,k>/ 2 

j ,k “ ( u i+*4, j+^.k + u i+*i,j-H,k + u i-k,j+H,k + u i-4, j -H,k)/ 4 
. . . etc. 

When extended to the other flux components and to averaging over cube 
surfaces in three dimensions, the numerical potential equation is of the 
form: 


| ^^(phU) + + - 0 

! 

! 
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To find the flux quantities />hU, phV, and phW at the finite volume 
cell vertices (i.e. points a, b, c t and d for the two dimensional case), 
it is necessary to evaluate Equations (2) through (4). The derivatives 
in these expressions can be expanded by the same volume averaging 
approach used above, thus: 

*, - *««,<•> ye - Nr 6 e ( y> 

z e “ M »7r 6 e (2) 

with similar terms for the other transformation metrics. The above 
expressions, being centered at grid midpoints, will involve the values 
of the potential and grid position at grid points which are known from 
the previous potential solution and the grid geometry, respectively. 

As mentioned in the Background discussion, when solving transonic 
flows it is necessary to include in the solution algorithm some form of 
supersonic upstream dependence in order to account for both the physical 
nature of the flow and the viscous nature of shock waves, respectively. 
Caughey and Jameson introduced upwinding by the addition of terms into 
their potential numerical equation which are only non- zero when the flow 
is supersonic. Also, the finite volume technique used exhibits a ten- 
dency for uncoupling of the flow field solution between alternating grid 
points. As a result, additional terms are included in the numerical 
potential equation. The final numerical equation which is solved by 
FL030 when these terms have been included has the form: 

^ r S f (phU+P) + + M^S^CphV+R) 
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where P, Q, and R are the upwinding terms and Q^, Q n j- . Qff- and Qfr/f 
are the decoupling terms. 

A reader interested in an in depth description of the nature and 
derivation of the additional terms mentioned above may find this inf or- 
mation in Reference 14. 

Computational Grid Geometry 

As mentioned previously, the computational grid used by FL030 is a 
body fitted, non - orthogonal , curvilinear mesh constructed about a wing 
fuselage combination. The number of grid points composing the computa- 
tional domain is typically 40 x 6 x 8 , 80 x 12 x 16, or 160 x 24 x 32 
for the number of £ , r>, and f points in the coarse, medium, and fine 
grids, respectively. The grid is conformally mapped to the wing and 
fuselage surfaces as can be seen from the plot of surface grid lines 
shown in Figure 3 . 

The grid is formed around spanwise airfoil sections in a similar 
manner in which "C" grids are mapped to airfoils in two-dimensional ana- 
lysis. This mapping for the two-dimensional case unwraps the physical 
plane from around the airfoil and from along a slit extending from the 
airfoil trailing edge to the downstream boundary. The wing surface and 
slit are then a line of constant t] as are the upper, lower and upstream 
boundaries. The upper and lower half outflow boundaries are lines of 
constant £. 

The FL030 grid differs from the "C" grids due to the need to map the 
grid to the fuselage surface in addition to the wing surface. As a 
result, rather than allowing the constant f planes to extend to the 
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upper and lower boundaries, these surfaces are wrapped around the fuse- 
lage and a line extending forward from the fuselage nose. The line of 
constant r, which forms the upper and lower surface boundaries in two- 
dimensions now becomes the upper and lower half plane symmetric bound- 
aries for the three-dimensional grid. Figure 4 shows a plot of the con- 
stant f surface which coincides with the wing tip airfoil section. Note 
that for this grid geometry, the far field side boundaries are all 
formed by the maximum r, surface; and the upstream and downstream bound- 

aries are the same as for the "C" grid case. 

In addition to the grid surfaces which form the physical extents of 
the computational grid, additional "ghost" coordinate surfaces are auto- 
matically generated below the wing and fuselage surfaces and beyond tne 
symmetric boundary. These "ghost" surfaces are necessary for the formu- 
lation of both the finite -volume numerical flow equations and the flow 
tangency boundary conditions upon these boundaries. The grid points 
composing the "ghost" surfaces are calculated from linear extrapolations 
of the computation mesh grid points immediately inside the physical 
domain. 

Boundary Conditions 

The original FL030 code included the following physical boundaries: 
solid surfaces, far field boundaries, symmetric plane boundary, and 
• trailing edge slit boundary. Note that when viscous effects are 

included by adding on to the wing surface a displacement thickness, the 
new displaced surface acts like a solid surface and has the same bound- 
ary condition as for the inviscid wing surface. However, as will be 
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discussed in the next section, inclusion of viscous wake effects will 
influence the boundary conditions applied across the trailing edge slit. 

Since the governing potential equations are written in terms 
of perturbations from free-stream conditions, the subsonic, far- field 
requirement that the flow return to the free-stream velocity and direc 
tion are satisfied by setting the perturbation potential equal to zero 
on the side and upstream boundaries. These boundaries are formed by the 
maximum f grid surface and part of the minimum i? surface as mentioned 

previously in the grid geometry section. 

For a purely subsonic and inviscid flow, the zero perturbation 
potential boundary condition would also exist on the downstream outflow 
boundary due to the isentropic (reversible) nature of the flow field. 
When solving transonic flows with shock waves or when viscous effects on 
the wing surface are being included, the flow field processes are irrev- 
ersible and a return to free-stream conditions will not be expected. The 
downstream boundary condition allows for these effects by simply utili- 
ing a "zero" order extrapolation of the potential (constant potential 
assumption) to the outflow boundaries. 

A flow tangency condition is applied along both the wing and fuse- 
lage solid surfaces by setting the normal contravariant component of the 
velocity vector to zero on the surfaces. This condition provides an 
equation which when approximated by a finite -difference expansion about 
the surface grid points can be used to set a value for the perturbation 
potential on the "ghost" grid points below each surface. Note that this 
finite -difference boundary condition is not exactly similar in form to 
the finite-volume solution algorithm of the governing equations. As a 
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result, it is possible to impose flow tangency using the finit- 
difference technique yet still have some normal surface velocity when 
performing the finite-volume calculations. Since it is essential to 
have the most exact boundary condition imposed upon the wing surface in 
order to generate the most accurate solutions, a second condition exists 
upon the wing surface. This additional condition involves reflecting 
the flux quantities calculated by the flow solver for the cell centers 
directly above the wing surface to the "ghost" cell centers beneath. 

The reflected normal fluxes then cancel each other out in the residual 
expression and a net zero flow is obtained through the surface. 

Since FL030 only solves the flow field for a half body configura- 
tion, only cases of no sideslip can be analyzed. Thus, the symmetric 
boundary has no flux through it and a flow tangency condition exists on 
the half plane. This condition is imposed numerically in the same fash- 
ion as for the wing surface with both the weak finite -difference and the 
strong finite -volume methods being used. 

For an inviscid calculation, the trailing edge slit boundary is not 
an actual limit to the physical domain as the other boundaries are but 
is simply an artificial boundary created by unwrapping the physical 
plane into the computational domain. The only conditions which need to 
be imposed at the slit is that the flow velocities, and thus pressure, 
be continuous across the cut. The flow potential, however, will have a 
discontinuous jump across the wake which is proportional to the sec- 
tional wing lift coefficient. 


/ 00 
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Viscous Boundary Layer and Wake Effects 

The program routines for the calculation of the viscous boundary 
layer and wake have been integrated into FL030 as part of an investiga- 
tion by'streett 17 into more accurate methods of correcting potential 
flow solutions for viscous effects. These routines model the three 
dimensional laminar and turbulent surface boundary layer with the tran- 
sition location specified by the user. The turbulent, viscous wake 
trailing the wing surface is modeled as having both finite thickness and 
curvature that vary with downstream location. The procedure for obtain- 
ing a converged viscous solution is to alternate between solving the 
potential flow and boundary layer equations until both solutions are 

converged. 

As mentioned previously, the effect of the surface boundary layer on 
the potential flow is accounted for by assuming a displacement thick- 
ness. These displacements are calculated by the boundary layer routines 
at the points which are originally input by the user for the wing geom- 
etry. When regenerating the computational grid after a viscous calcula- 
tion, the displacements are added normal to the original surface and the 
grid formed around this displaced surface. The boundary conditions on 
the displacement surface is the same as for a physical surface, i.e. 
flow tangency. 

The wake thickness correction is similar to that for the surface 
displacement thickness and is also included in generating the computa- 
tional grid in order to produce a trailing edge slit which has a non- 
zero width. However, the boundary conditions applied across the finite 
different from the inviscid potential flow case due to the pos- 

(01 


slit are 
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sibility of a pressure jump across the wake and, thus, some wake curva- 
ture in the downstream flow. The presence of wake curvature is imposed 
on opposing sides of the wake cut, as a difference in potentials, which 
decays to zero or freestream conditions far enough downstream. 

Since, the present investigation deals with the design of the wing 
surface itself, alternate boundary conditions on the design surface will 
be required. In the existing analysis code, only the general potential 
flow solution method, the wing/displacement boundary conditions, and the 
application of the boundary layer displacement thickness on the surface 
of the wing are of immediate concern in this thesis. 


/OX 
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INVERSE WING DESIGN METHODS 

As stated previously, a direct -inverse approach to wing design was 
selected for incorporation into the TAVFIVE code. The direct- inverse 
method derives its name from the division of the design wing surface 
into a fixed geometry leading edge region, where flow tangency boundary 
conditions are imposed, and an aft, variable geometry section where 
pressure boundary conditions are enforced. The pressure boundary where 
the user specified pressure distributions are imposed does not extend 
forward to the leading edge due to difficulties of enforcing this type 
boundary condition near the beginning of an airfoil section. This 
restriction on the size of the pressure specification region does not 
seriously reduce the versatility of the design method since the leading 
edge regions for most airfoils are similar, and it is relatively easy to 
select a leading edge geometry which will produce the desired Mach num- 
ber or pressure values at the beginning of the inverse region. In addi- 
tion, specific leading edge shapes may be required due to other design 
constraints such as the necessity to house a leading edge flap or slot 
system. Finally, as will be shown later, the application of the inverse 
approach only to back portions of the wing aft of the leading edge means 
that it is used in regions where the pressures are primarily dominated 
by the chordwise flow velocity; and, thus, the direct- inverse approach 
allows some useful simplifications when formulating the equations. 

The following sections will provide details concerning the usage of 
the pressure distribution as a boundary condition in the inverse wing 
regions and how a new surface which would produce these pressures can be 
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calculated. An additional feature which helps to restrict the inversely 
designed wings to practical geometries by imposing a desired trailing 
edge thickness is also described. The final section details how the 
inverse design procedures are integrated into the program as a whole and 
presents a logical design strategy for use in designing actual wings. 
Part of this last section also discusses the code's data input structure 
and design options available. 


Pressure Boundary Condition 

In the inverse design regions on the wing, a pressure boundary con- 
dition will be specified rather than the flow tangency condition used -in 
analysis zones. In formulating this boundary condition it is necessary 
to relate the user specified pressure coefficient, Cp, to the current 
perturbation potentials at inverse design grid points. Consider the 
full potential equation for the pressure coefficient: 


where 


r ^ 


C P " 2 i 


7-1 


1 + 2 • 2 ^ 
q® 
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- 1 




+ w* 


If it is assumed that the pressure coefficient is primarily a function 
of the chordwise component of the velocity, u, and only slightly 
affected by the vertical and spanwise components of velocity, v and w, 
then the latter two velocities can be time lagged in the boundary condi- 
tion without introducing any solution Instabilities . This assumption is 
true everywhere except near the leading edge; but since the inverse 
design boundaries have already been restricted to regions well behind 
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the leading edge, the simplification is justified. The value of the 
local velocity, u, can then be calculated from the above expression in 
terms of the desired pressure coefficient and the current values for the 
vertical and spanvise velocities. In addition, the velocity u can also 
be calculated from the perturbation potentials using the relations of 
Eq. (2). Defining to be the elements of the inverse transpose of 

the Jacobian matrix, H, the two equations for u yield: 


1 

J ll^£ + J 12^»? + J 13^f “ *" 



cos (o ) (5) 


Since the spanwise and vertical flow velocities have already been 
assumed to be constant in the boundary condition relation, it is consis- 
tent to make the same approximation in the above expression with respect 
to the spanwise and vertical derivative terms, ^ and <f>^ . This 
assumption is similar to the previous one, and it will lead to an expli- 
cit expression for the potential at one point. 

Two different finite difference approximations of Eq. (5) have been 
considered. The first method involves expanding the derivatives of the 
potential function as central differences about a grid point in the £ 
(chordwise) and f (spanwise) directions and a three point backwards 
finite difference in the * (vertical) direction. The reason for the 
special treatment on the tj derivative is due to the value of the poten- 
tial not being defined at the "ghost" points in the inverse regions. 

The numerical equation with the above expansions is. 
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n+1 n 

,k ' ,k)/ 2 

n n n 

+ ,k * 4 *i-l,j-l,k + *i-l,j-2,k>/ 2 

+ J 13<^i-1 . j ,k+l * ,k-l>/ 2 “ F ( c Pi-l,k> 

Here, the superscripts n and n+1 refer to current values of the 
potential and the new values of the potential being imposed by the 
boundary condition, respectively. Also, the term F(Cp£.^ is the 
right hand side of Eq. (6) evaluated using the pressure coefficient spe- 
cified at point i-l,k. Figure 5 shows the points involved in the above 
expression and the point at which the pressure is being specified. This 
equation tan then be solved for the potential at grid point i+l,j,k 
which yields the equation: 

n+1 _1_ n 

" J n J ll*i-2,j,k 

n n n 

- <*120*1-1, j.k * ^i-l.j-l.k + *i-l,j-2.k) 

- J 13 W 1 -IJ ,k+l - *i-l,j,k-l) + 2F ( c Pi-l.k)J 

The potential values at n+1 in the direct region are known initially 
since they do not change when the inverse boundary condition is applied; 
i.e. - <^ n . All the potentials on the. inverse boundary can then be 

calculated and, since the spanwise and vertical derivatives are small, 
W 1H primarily be functions of the pressure coefficient at grid point 
i-1 and the value of the potential at grid point i-2. 

This method of approximation has the possibility of producing two 
independent solutions for the alternating set of even and odd grid 
points which is a problem commonly seen when using simple central dif- 
ference approximations. The reason that this solution uncoupling may 
exist is due to the dependence of each potential on the value two grid 
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Figure 5. Point Dependance for Grid Point Specification Method 
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points ahead or behind rather than the one immediately adjacent. To 
preclude the possibility of this problem occurring, another finite dif- 
ference method has been considered. 

The second finite difference approximation involves expanding the 
derivatives of the potential about the mid-point i+4, j ,k. The £ deriva- 
tive for this case is determined by a central difference involving the 
preceding and following grid point values. The rj and f derivatives can 
be found at the mid-point by averaging the derivatives from the preced- 
ing and following grid points found be the three point backwards and 
central difference approximations as in the previous method. Figure 6 
shows the point dependence and pressure specification point for this 
method. The resulting numerical expression obtained with these finite 
approximations is : 

n+1 n 

t , ,, n+ l n n n 

+ J 12 + *i-l,j,k> ' 4 <*i,j-l,k + *i-l,j-l,k> 

n n 

+ <*i,j-2,k + *i-l,j-2,k / 4 

n n n n 

+ Ji3(0 it j ,k+l + *i-l,j,k+l * *ij,k-l * *i-l,J,k-l>/ 4 

" F < c Pi-h,k> 

from which the potential value at i+l t j ,k can be found as: 
n+1 • 

1 

J ll + 3J 12/ 4 

* J 12 3 ^i-l,j,k - 4 <*i,j-i,k + *i-l,j-l,k> 
n n 

+ ^i,j-2,k + ^i-1 , j -2 ,k / 4 
J 13 , j ,k+l + ^i-1, j ,k+l * ^i, j ,k-l * ^i-1 , j ,k-l)/ 4 _ 

+ F (Cpi..4 ik ) 
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Figure 6. Point Dependance for Mid-point Specification Method 




The only concern with using this mid-point specification scheme is 
that the current method of calculating the pressure data output from 
FL030 uses the first difference scheme for the streamwise derivative. 
This difference could potentially allow a pressure to be specified cor- 
rectly but still have a significantly different value output from FL030 
due to the inconsistent calculation methods. However, as shown on Figure 
7, where the pressures calculated for a typical flow solution are com- 
pared for the two different calculation techniques, this possible error 
has not been significant in practice. 


Surface Calculations 

As the inverse boundary conditions drive the flow field to a con- 
verged solution, it is necessary to periodically calculate the location 
of the new displacement surface and to -regenerate the computational grid 
about this new geometry so that the pressure boundary surface will cor- 
respond to the physical boundary surface. Each new surface can be found 
relative to. the previous surface from an integration of the wing surface 
slopes. However, the surface slopes must first be calculated from the 
current flow field solution using the flow tangency boundary condition 
which in Cartesian coordinates is: 

U T x VF - 0 

where U is the Cartesian velocity vector and VF is the gradient of the 
surface function F. 


Since all of the numerical calculations are performed in the curvil- 



inear computational space, it is necessary to express the boundary con- 
ditions in the same manner. Defining V'F to equal the gradient of F 
with respect to the curvilinear coordinates and V to be the contravar- 
iant velocity vector, the following relations can be obtained from Equa 
tions (2) and (3): 

U - H x V VF - (H* 1 ) 7 x V’F 

From these the transformed boundary condition can be obtained as: 

(H x V) T x [(H _1 ) t x V'F] - V T x [H' 1 x H] T x V'F - 0 

or simply 

. V x V'F - 0 

As can be seen, the tangency boundary condition written in terms of 
the curvilinear coordinate system using contravariant velocities is a 
direct analog to the same condition expressed in physical space. 

A more useful expression can be obtained by expanding the above 
equation to: 



This expression can be solved for the new chordwise airfoil slopes, 
if the current values of the spanwise slope, drj/d ^ , are used. 
Since the wing surface is represented in the computational grid as a 
plane of constant rj , the current slopes on the wing surface equal zero 
and a simplified flow tangency condition results: 


(ft). 
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This flow tangency condition is only exactly true when the veloci- 


ties are calculated at the new physical surface boundary. Using these 
velocities is difficult in practice, however, since the new surface 
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being calculated will typically lie at locations between grid points, 
where velocities are difficult to calculate; or even at locations out- 
side the flow field altogether, as when the new surface falls below the 
old surface. Thus, for simplification, the contravariant velocities 
used in calculating the surface slopes are the velocities at the inverse 
pressure boundary; i.e. the current grid boundary which is also the pre- 
vious physical surface boundary. This approximation is justified by the 
fact that for the converged surface solution the pressure and physical 
surface coincide and the equation becomes exact. In addition, no sur- 
face convergence problems resulting from this approximation have been 
observed, even for cases where the pressure and physical boundaries ini- 
tially differed by a significant amount. 

A second factor affecting the slope calculations is the accuracy of 
the contravariant velocity calculations. A first approach to calculat- 
ing the U and V velocity components was simply to apply the relations of 
Eq. (2). The derivatives of the potential function were evaluated by 
central difference approximations in the ( and f directions and a three 
point backwards difference in the r, direction. The later was required 
because, as has already been mentioned in the pressure boundary section 
of this thesis, the potential function was not defined at ghost points 
in the inverse regions. As was expected from the previous experience of 
Weed et al.^, this approach did not yield accurate values for the V 
velocity component. As can be seen from the data plotted in Figure 8, 
the surface displacements, labeled "Finite Difference Approach," based 
upon velocities calculated using the above method for. a converged analy- 
sis solution of a wing section are significantly different from zero. 
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Figure 8. Comparison of Accuracy cf Finite-Difference Method 

and Residual Method 
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The later is, of course, the proper displacement expected for a con- 
verged solution. 

The other data plotted on Figure 8, labeled "Residual Approach" , are 
the result of using an alternate approach which is an extension of the 
formulation of Weed et al. 9 . This second method involves solving the 
residual expression from the potential flow solver for the surface velo- 
city ratio, V/U, under the assumption that the residual is zero. For 
this problem, the residual expression from FL030 can be written in 
finite volume form as: 

Mj7f 5£(phU) + + ^5 r (phV) + (other terms) - 0 

The "other terms" in the above expression involve the grid point 
coupling and upwind dependence terms of the formulation and will be 
assumed to be constants in the following development. 

The desired velocity ratio, V/U, can also be written in this finite 
volume form as : 

V _ £hV _ ^ nr (phV) 

u phu my 

By simple manipulations, this ratio can be obtained from the residual 
expression as: 

H£ V f(phV) 2^^^(phV)^, 1 + MqfSfCphU) + **£ ^(phW) + (other terms) 

'2jx^ f (P hu >. 

where the subscript r? - 1 refers to the values at grid cell centers above 
the wing surface . 

In order to use Eq. (6) to find the desired surface velocity ratio, 
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it is necessary to know the U and V velocity components at the "ghost” 
cell centers below the wing surface. These values can be obtained in a 
manner consistent with FL030 by specifying the "ghost" cell values to 
equal the values at corresponding points immediately above the wing sur- 
face. Eq. (6) then explicitly defines the velocity ratio at the bound- 
ary grid points and, as seen from Figure 8, is very accurate. 

It should be noted that an interesting program simplification has 
been achieved in performing this calculation. In the discussion of the 
analysis surface boundary condition, it was stated that the "ghost" cell 
velocities were obtained by reflection about the wing surface. This 
procedure amounts to setting the chordwise and spanwise velocities, U 
and W, at "ghost" points equal to those above the wing surface and set- 
ting the vertical velocity component, V, at the "ghost" point equal in 
magnitude but in the opposite direction to the V velocity above the wing 
surface. If the residual were calculated in this manner at the surface 

i 

grid points in the inverse regions, then the surface slopes can then be 
found simply by: 

(iJL\ - nr (phV) - /Residual) (7 ) 

/wing 2^ r (phU) 

Thus, a converged surface solution where dr\/b£ tends to zero corresponds 
to a converged finite -volume solution on the surface where the residual 
tends to zero . 

With the contravariant velocities known, an integration of Eq. (7) 
through the inverse design region from the leading edge to the trailing 
edge will yield a set of surface displacements for the new wing surface 
relative to the previous one. These changes will be expressed as 


\\\c 



changes in the computational coordinate rj, which can be converted to 
surface displacements in the physical plane via the transformation: 

Ax - Xjj Atj Ay - y^ A*? 

However, since the computational grid lines are orthogonal at the 
surface, the normal displacement, 6(x), is: 

6(x) - ((Xjj hr]) 2 + (y n hr ] ) 2 ) 1/2 

These displacements will be defined at the computational grid points 
in the inverse regions. To obtain the corresponding displacements at 
the original geometrical locations specified in the program input data, 
a linear interpolation of the above data is performed. The reasons that 
the displacement surfaces are needed at the original geometry points is 
two fold. First, as mentioned previously in the discussion of viscous 
interaction, the boundary layer displacement thickness, which is numeri- 
cally analogous to the inverse displacement thickness , is defined at the 
original input stations: and finding the inverse displacement thickness 
at the same locations allows the use of the same routines for adding the 
boundary layer and inverse design displacements to the original geom- 
etry. Second, finding the displacements at the original geometry sta- 
tions permits the calculation of the new wing airfoil sections at the 
same semispan locations. 

Trailing Edge Closure 

The procedures outlined above will compute a wing surface corre- 
sponding to a given, fixed, leading edge geometry and to a desired set 
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of pressure distributions in the inverse regions. The above procedures 
do not, however, guarantee that this wing geometry will be practical. 

In particular, past experience^has shown that inverse surface calcula- 
tions may yield airfoil sections which have either excessively blunt 
trailing edges or which, at least numerically, have the upper and lower 
surfaces crossed at the trailing edge ("fish tailed"). The former case 
is undesirable due to aerodynamic considerations, while the latter is 
physically impossible and may produce unpredictable problems in the grid 
generation or flow calculation portions of FL030. 

Since for any specified pressure distribution the corresponding 
wing surface will be controlled by the leading edge geometry, which 
serves as an initial spatial boundary condition for the inverse region, 
the problem of assuring trailing edge closure can be viewed as the 
proper selection of a leading edge shape. A procedure for systemati- 
cally modifying the leading edge region in order to achieve some desired 
trailing edge thickness is called relofting. Such a relofting procedure 
has been incorporated into the present design process in order to both 
prevent the problems of trailing edge crossover and to allow the user 
the option of specifying a trailing edge thickness as an additional 
design variable. This design feature should be very useful in practical 
applications since it automates the Iterative selection of a leading 
edge shape which would otherwise have to be performed by the user. 

The present method of surface relofting is a simple linear rotation 
scheme. This method can be visualized with the help of Figure 9. The 
dashed line indicates the original leading edge geometry and a hypothet- 
ical new surface shape which has been calculated for the inverse design 
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regions. Without modification, this new surface has a trailing edge 
thickness of A. If a thickness of A t was specified by the user, then 
the surface would have to be relofted or changed. In the present 
scheme, in order to obtain the desired thickness, a displacement thick- 
ness, £ r , is added to the current design surface. This thickness has a 
distribution from the leading to the trailing edge and is determined by 
the formula: 

£ r (x) - (A t - A) (x/c) 

where c is the chord length of the local airfoil section. The total 
displacement for a surface update is then the sum of the two displace- 
ments, 6 and £ r . When both the upper and lower surfaces are designed 
simultaneously , the displacement* magnitudes determined by relofting are 
divided between the two surfaces so that half is added to the lower sur- 
face and half to the upper surface. 

It should be noted that the above procedure for relofting the lead- 
ing edge region is a practical method, and it is not based upon any 
known dependence between the leading edge shape and the trailing edge 
thickness. Experience has shown, however, that this method does indeed 
produce the appropriate effect at the trailing edge; and with iteration, 
a converged geometry solution with the desired trailing edge thickness 
can be obtained. 

Design Strategy 

All of the above procedures are part of the overall wing design pro- 
gram structure. The addition of the necessary program statements and 


subroutines to the TAWFIVE code has been made in an unobtrusive manner 
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so that the original TAWFIVE program logic and input formats are pre- 
served. This section will discuss in detail how these inverse design 
methods interface with the analysis portions of the code and how 

inverse design of a wing would be performed. 

Consider first how the TAWFIVE program is used when performing an 
analysis of a specified wing geometry. The wing and fuselage geometry 
data are input and the desired computational grid is subsequently gener 
ated. Typically, there is a choice of three different computational 
grids which may be thought of as coarse, medium, and fine. The usual 


procedure is to start a wing analysis on the coarse grid and then to 
halve the grid spacing to the medium size after some specified number of 
iterations of the potential flow solver have been performed. After 
additional potential flow iterations on the medium grid, the grid spac- 
ing is again halved to the fine grid size; and. for purely inviscid ana- 
lysis. iterations are repeated until a converged potential flow solution 

is obtained. 

For an analysis with viscous interaction, it is sufficient y 

partially converge the potential solution on the fine grid so that rea 
sonable wing surface pressures can be calculated; and the boundary layer 
routines are then called in order to calculate an initial displacement 
thickness. With the displacement surface known, the procedure generates 
a new grid around the displaced surface; and the iterative and grid hal- 
ving procedure is repeated until a converged solution for both the 
potential flow and the boundary layer is obtained. Note that as the 
solutions converge, it is usually not necessary to return to the coarse 
grid after each displacement surface update; but solutions can continue 
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on the fine grid starting with the most recent potential solution. 

As shown from the flowchart in Figure 10, an inverse wing design 
will proceed in a very similar manner to that of a viscous wing analy- 
sis. After the initial wing geometry is input and the computational 
grid has been generated, the inverse pressure distribution data and 
design control parameters are inputted and the pressures to be specified 
at the grid mid-points are calculated. Before each potential flow iter- 
ation is performed, a call is now made to the pressure boundary subrou- 
tine where the potentials on the wing surface in the inverse regions are 
specified. The potential flow solver then proceeds as normal except 
that the changes of the potential values computed on the inverse wing 
surface are set to zero. After all the desired grid halvings are per- 
formed and a semi -converged solution is obtained, the surface update 
routines are called. These routines compute both a new surface geometry 
and a set of inverse displacement thicknesses. The new geometry is out- 
put to a print file in the same format as the input geometry file and is 
not saved in memory. The displacement thicknesses , however, are 
retained and are used in computing a new computational grid in the same 
manner as the viscous displacement thicknesses are used. 

When performing an inverse design with viscous interaction effects 
included, the procedure is exactly the same as the inviscid design 
procedure with the viscous and inverse displacement thicknesses being 
used in the direct and inverse regions of the wing, respectively. It 
has been observed, however, that it is best to obtain a semi -converged 
boundary layer solution (at least one iteration) before beginning the 
inverse design. The reason for needing an initial boundary layer solu- 
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Figure 10 . Flowchart of Inverse Design Procedure 
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tion is due to an occasional over prediction of the boundary layer 
thickness at the trailing edge by the first boundary layer calculation. 
These excessive, erroneous thicknesses may cause the physical surface 
beneath the displacement surface to appear to be fish tailed. If forced 
trailing edge thickness is being imposed, the result is an unnecessary 
surface relofting which must be undone by additional surface updates, 
thus slowing down the solution convergence. 

Design Input and Control 

The inverse pressure file is oriented around individual span sta- 
tions in the same manner as the wing geometry input file is organized 
about sectional airfoil geometries. As such, pressures are input as 
chordwise distributions of the pressure coefficient for each span sta- 
tion with as many span stations as necessary being specified in order to 
fully define the desired wing pressure distribution. Every set of two 
spanwise stations may be thought of as a design panel on the wing sur- 
face, with each panel being either an upper surface design, lower sur- 
face design, or both. Since input pressure stations and computational 
grid span stations will not necessarily correspond to one another, once 
a design panel is specified, the program searches to find which, if any, 
grid stations are positioned within the inverse region. Pressures and 
the other design variables are then linearly interpolated from the sur- 
rounding pressure stations to obtain the values at the grid station. 
These other variables are the location of the direct- inverse junction 
near the leading edge and the trailing edge thickness desired for trail- 
ing edge closure. Contiguous design panels which involve design of the 
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same surfaces (i.e. upper, lower, or both surfaces) can be formed by 
inputting a single additional pressure station which will define the 
panel between the last input station and the new. Figure 11 shows the 
input pressure sections and the resulting computational inverse regions 
for a design case which involves both a continuous design region defined 
by multiple pressure input stations and a separate design section near 
the wing tip. The computational grid lines shown are typical for medium 
grid spacing. Note that the design region between the first two pressure 
stations does not contain any computational grid lines and thus will 
have no affect on the final wing design. This example indicates that 
the user should be familiar with where the computational grid will be 
formed before arbitrarily selecting input stations. 

Two control parameters have been added to the TAWFIVE iteration con- 
trol input file. The first controls whether or not the inverse boundary 
conditions will be imposed for a given set of potential flow iterations 
and, if the conditions will be imposed, at which iteration count they 
begin. This latter option has been included because it is foreseen that 
for geometry cases with a difficult convergence trend, it may be neces- 
sary to obtain a semi -converged solution before attempting to enforce 
the design conditions. All of the design cases which will be discussed 
in the results section of this thesis have been obtained by enforcing 
the design boundary conditions starting with the first iteration and no 
convergence difficulties have been observed. 

The second control variable is used by the surface updating routine 
to determine what type of relofting control is desired. Since relofting 
will always occur it the current wing surface will result in a computa- 
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tional grid which crosses itself prior to the trailing edge, this vari 
able is used only to determine if the trailing edge thickness is to be 
"forced" to a user specified value. Note that relofting will automati- 
cally occur if the grid itself crosses, but it will not take place if 
only the airfoil crosses itself. The distinction is due to the fact 
that when the viscous boundary layer is added to the problem, it is pos- 
sible that the airfoil may be fish tailed while the displaced surface 
not. In these cases a fish tailed airfoil may actually be designed if 
the user does not force trailing edge closure. 
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RESULTS 


In this section results from five different test cases for both sub- 
critical and supercritical conditions will be presented. These cases 
are not intended to be definitive or even representative of practical 
designs but have been selected as examples of the capabilities of the 
present inverse design technique. The results shown were obtained on a 
medium grid having 81 streamwise, 13 vertical, and 19 spanwise points 
with 11 spanwise stations and 53 points on the wing at each station; and 
in all cases the maximum change in the reduced potential was reduced at 
least three orders of magnitude. Thus, the results do not represent 
ultimate convergence but should be representative of "engineering accu- 
racy" . 

The planform selected for the test cases was the Lockheed Wing A 
wing-body. The wing for this configuration has a quarter chord sweep of 
25 deg., a linear twist distribution ranging from 2.28 deg. at the wing 
body junction to -2.04 deg. at the wing tip, an aspect ratio of eight, 
and a taper ratio of 0.4. The last two values are based upon the wing 
without fuselage. However, instead of the supercritical sections nor- 
mally associated with Wing A, the initial airfoil sections at each span 
station were assumed to be NACA 0012 airfoils. 

The target pressure distributions used in the design regions were 
selected to yield airfoil shapes thicker in the aft portions of each 
section; and, at supercritical conditions, to yield on the upper surface 
weaker and more forward shock waves than those which would normally 
occur on a NACA 0012 section. On the lower surface, the target pressure 
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distributions were selected to have either a favorable pressure gradient 
or fairly constant pressure plateau over much of the lower surface. 

All subcritical cases were for a freestream Mach number of 0.7 and 
an angle of attack of two degrees. In each case, the pressure distribu- 
tion was specified in the design regions from the 15% local chord loca- 
tion to the trailing edge and used as the boundary condition in these 
inverse regions starting with the first iteration. Normally, two 
hundred SLOR iterations were executed prior to the first design surface 
update calculation; and subsequently, surface updates were computed 
every fifty cycles. Usually, the solution was considered converged and 

terminated after 450 total iterations. 

Supercritical cases followed a similar procedure except that the 
freestream Mach number was 0.8. Again the angle of attack was two 
degrees. However, for these cases three hundred iterations were per- 
formed prior to the first surface update calculation in order to better 
resolve the leading edge pressure distribution in design regions. 

Because of the upstream dependance of the flowfield, particularly for 
the supercritical cases, it was determined to be essential to obtain a 
good computational solution in the leading edge region before any sur- 
face updates. Otherwise the initial surface changes were so drastic 
that a large number of additional surface calculations, and accompanying 
iterations, were necessary in order to achieve convergence. 

Finally, for those cases where trailing edge closure was specified 
by the input, forced relofting was not performed until the second sur- 
face update. This approach was used because the first surface update 
usually involved large changes in the surface shape, and it was believed 
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that attempting to force closure at the same time might lead to conver- 
gence difficulties. 

Test Case A 

As shown on Figure 12, the objective of Case A was to modify only 
the upper surface between 45% and 85% semi-span. As indicated above, 
the input pressure distribution for the design region corresponded to 
that of a wing composed of airfoil sections which were thicker than a 
NACA 0012 in the aft portion of each section; and these pressures were 
previously obtained with a corresponding analysis computation. Thus, 
since this case also required trailing edge closure, Case A was a test 
of the ability of the method to reproduce the airfoil sections of a 
known wing. Both subcritical, designated Case Al , and supercritical, 
designated Case A2 , solutions were obtained. 

The resultant designed airfoil sections for the case having a sub- 
critical freestream are shown on Figure 13. As can be seen, the 
designed sections are considerably different than the original NACA 0012 
airfoils; and they are in reasonable agreement, even on the expanded 
scale, with the target sections. However, there are some slight discre- 
pancies at the boundary stations at 50% and 80% semi- span. It is 
believed that these are due to a combination of terminating the computa- 
tion prior to ultimate convergence and to the significant variation in 
spanwise slope near the trailing edge resulting from the change between 
the NACA 0012 sections in the analysis zones to the designed airfoils in 
the inverse regions. Nevertheless, it is believed that the agreement 
between the designed surfaces and the target surfaces is adequate. 
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The true test, however, of an inverse wing design method is not its 
ability to reproduce "known" airfoil sections but rather a comparison 
between the target pressure distributions used to design the wing and 
those computed by an analysis of the designed wing. Figure 14 presents 
such a comparison for subcritical Case Al; and, as can be seen, the ana* 
lysis results for the designed wing (labeled "designed surface pres- 
sures") are in excellent agreement with the target pressures as are the 
local lift coefficients. 

Figure 15 and 16 show similar section profiles and pressure distri- 
bution for Case A2 at supercritical conditions. Again the agreement 
between the designed surfaces and the target surfaces and the pressures 
from an analysis of the designed wing and the target pressures are 
excellent. It is believed that Figures 13-16 demonstrate that the cur- 
rent method can be used to modify the design of the upper surface of a 
wing mounted on a body. 

Test Case B 

This case, which is depicted on Figure 17, was created to test the 
ability of the method to design both upper and lower surfaces. Subcrit- 
ical (Case Bl) and supercritical (Case B2) results are shown on Figures 
18-21. As in the previous case, trailing edge closure was required; and 
as a result the designed surface shapes have the same character as those 
for Case A in that there is good agreement at the inner stations but 
slight discrepancies between the designed surfaces and the target sec- 
tions at the boundary stations. However, as shown on Figures 19 and 21, 
there is still excellent agreement between the pressures computed by an 
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analysis of the designed wings and the desired target pressures used in 
the inverse design. Thus, it can be concluded that the method can be 
used to modify the design of the upper and lower surfaces of a wing 
mounted on a body. 

Test Case C 

The inverse design regions for Case C, which was an attempt to 
design both upper and lower surfaces on two noncontiguous regions of the 
wing at supercritical conditions, are shown on Figure 22; and a compari- 
son between the initial pressure distribution associated with NACA 0012 
sections and the target pressures is portrayed on Figure 23. As can be 
seen, the target pressure distribution essentially eliminates at inboard 
stations the upper surface shock wave present on the original wing; and 
at outboard stations it weakens the shock and moves it forward. In 
addition, significant changes in the lower surface pressure gradients 
are evident. Also shown on Figure 23 are the pressures computed by the 
program at the end of the inverse design procedure (denoted as "design 
pressures"). These pressures are in excellent agreement with the target, 
pressures, which indicates that the method is satisfying properly the 
desired inverse boundary conditions. 

The corresponding designed airfoil sections for this case are shown 
on Figure 24. Even on the expanded scale, the agreement between the 
designed and target surfaces is excellent at all design stations. How- 
ever, trailing edge closure was not enforced for this case; and there is 
at the boundary stations some departure between the designed, surfaces 
and the target surfaces near the trailing edge. Again it is believed 

is-o 


71 



Figure 22. Inverse Design Regions for Case C 
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Figure 23. Comparison of Initial Pressures with Target Values (Case C) 
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that this slight difference is a ramification of the change in spanwise 
slopes near the trailing edge between the direct and inverse regions . 

In any event, the pressure distributions resulting from an analysis 
of the designed surfaces shown in Figure 24 are in excellent agreement 
with the target pressures, as can be seen on Figure 25. In addition, 
the section lift coefficients at the various design stations are in very 
good agreement with the target coefficients. Based upon these results 
it is believed that the present method can adequately design/modify non- 
adjacent regions of a wing in transonic flow. 

Test Case D 

As shown On Figure 26, Case D involved the inverse design of the 
entire wing on both the upper and lower surfaces. In addition, as 
depicted on Figure 27, the initial twist distribution was constant from 
root to 40% semi -span followed by a linear distribution between 40% and 
the wing tip; and the inverse pressure distribution was selected to cor- 
respond to an approximately linear twist distribution between the root 
and the tip. Thus, this case was a test of both the ability of the 
method to design an entire wing and to modify the twist distribution. 
Obviously, since the twist had to be permitted to vary, trailing edge 
closure was not required. Also, the results shown are for supercritical 
conditions . 

As can be seen on Figure 27, the twist distribution resulting from 
the design calculation, while considerably different than the initial 
distribution, is slightly different than the target distribution. This 
difference occurred for several reasons. First, in the current version 
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of the program the wing section at the root-body junction cannot be 
inversely designed. Thus, when designing the entire wing, the program 
automatically makes the root section nondimens ionally identical to that 
at the first span station; and the twist at the root and at the 10% 
semi-span station are identical. Second, the leading edge shapes in the 
direct region forward of 15% chord correspond to the initial shapes and 
are oriented by the initial twist distribution. Thus, they do not cor- 
respond to those associated with the target twist. Consequently, if the 
method correctly matches the input pressure distribution in the inverse 
region from 15% chord aft, it should yield slightly different pressures 
near the leading edge and a slightly different final twist distribution. 

Figure 28 compare the designed airfoil sections with the original 
surfaces. Due to the manner in which these plots were constructed, if 
the trailing edge of a designed surface is above that of the correspond- 
ing original surface, then that design station has a lower twist angle 
than the initial twist. As can be seen from Figures 27 and 28, the 
designed wing is considerably different than the original and has an 
almost linear twist distribution. 

As indicated above, the only way a design can be validated is to 
analyze the designed wing and compare the resultant pressures in the 
inverse regions with the target values. Figures 29 present such a com- 
parison for Case D, and it is apparent that the present direct- inverse 
method did design a wing having the appropriate pressures in the inverse 
regions aft of 15% chord. However, as should be expected, since the 
leading edge regions were different then those. corresponding to a true 
linear twist case, the pressure distributions in the leading edge 
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regions and the section lift coefficients were slightly different than 
those of the target case. (The target lift coefficients were obtained 
by an analysis of the target section shapes with a linear twist distri- 
bution.) It is believed that the results shown on Figure 29 demonstrate 
that the present method can be used to design an entire wing in super 
critical flow. 

Test Case E 

As a final test case, it was decided to design two non- adjacent 
upper surface regions simultaneously with a lower surface region which 
overlapped the upper zones. The location of these inverse design 
regions is shown on Figure 30. Likewise. Figure 31 compares the pres- 
sures associated with the initial wing sections shapes to the target 
pressures and to the pressures computed at the end of the design calcu- 
lation. It should be noted that this case is for supercritical condi- 
tion and trailing edge closure is not enforced. As can be seen, at sta- 
tions where only one surface is being designed (e.g. 20%, 40%, 50%, ana 
70%) the pressure distribution on the fixed surface also changes due to 
three dimensional effects from adjacent station which have been rede- 
signed. However, as depicted on Figure 32, only the design surfaces 
change form the original shape; and these surfaces are in reasonable 

agreement with the target profiles. 

Finally, Figure 33 compares analysis results obtained for the 
designed wing with the target pressures. Even for this complicated 
case, the agreement between the two distributions and between the actual 
and target lift coefficients is excellent. 











Figure 31. Comparison of Initial Pressures with Target Values (Case E) 
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CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK 

A direct- inverse wing design method has been successfully incorpor 
ated into the TAWFIVE transonic wing-body analysis computer code. The 
resultant code is capable of designing or modifying wings at both tran- 
sonic and subsonic conditions and includes the effects of wing-body 
interactions. A series of test cases have been presented which demon- 
strate the accuracy and versatility of this inverse method. 

Inclusion of viscous effects via the addition of the wing surface 
displacement thickness and wake thickness when performing wing design 
has been accomplished but not completely verified. Additional work will 
be required to run a sufficient sampling of test cases for evaluation of 
this design mode. The unique problems associated with viscous design 
and the effects of the various viscous correction models available m 
TAWFIVE would be the subject of a continuing research effort. 

The development and evaluation of alternate methods of surface 
'relofting are also topics for which continued research is suggested. The 
current method of relofting restricts the user to a family of leading 
edge geometries which can be constructed by the linear rotation of the 
initial shape. The option of using other relofting methods would extend 
the family of available shapes and add versatility to the design method. 
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abstract 


Verification, Optimization and Refinement of a Direct-Inverse Transonic 
Wing Design Method Including Weak Viscous Interaction. (August 1989) 
Robert R. Ratcliff, B.S., Texas A&M University 
Chair of Advisory Committee: Dr. Leland A. Carlson 


New developments in the direct-inverse wing design method in curvilinear co- 
ordinates are presented. A spanwise oscillation problem and proposed remedies are 
discussed. Test cases are presented which reveal the approximate limits'on wing as- 
pect ratio and leading edge sweep angle for a successful design, and winch show the 
significance of spanwise grid skewness, grid refinement, viscous interaction, the initial 
airfoil section and Mach number - pressure distribution compatibility on the final 
design. Furthermore, preliminary results are shown which indicate that it is feasible 


to successfully design a region of the wing which begins 


aft of the leading edge and 


terminates prior to the trailing edge. 



IV 


To Mom 


/I* 



ACKNOWLEDGMENTS 


The work presented in this paper was primarily supported by the National Aero- 
optics and Space Administration under Grant NAG-1-619 with Richard L. Campbell 
ol the Langley Research Center as technical monitor. The author expresses his appre- 
ciation to Dr. Carlson for his support, guidance, and patience and to Thomas Gaily 
for his assistance and helpiul suggestions. Much thanks also goes to the Computer 
Service Center at Texas AfcM University for their very generous computer support. 


I H 1 ! 



NOMENCLATURE 


A 

AR 

COS 

c 

Cl 

Cl 

C P 

°? 

cosh 

F 

a 

a i,j 

d 

f 

\h\ 

H 

h 

I,J,K 

J 

M 


influence coefficients used in compensation terms 
aspect ratio 
cosine 
local chord 

airfoil section lift coefficient 

wing lift coefficient 

pressure coefficient 

specific heat at constant pressure 

hyperbolic cosine 

Wing surface function in the physical domain 
speed of sound 

Fourier coefficients used in grid scheme 

the relative x distance from the sectional quarter chord point 
general function 

determinant of the inverse jacobian matrix 
inverse Jacobian transformation matrix 
enthalpy per unit mass 
grid locations in £, 7?, C directions 
Jacobian transformation matrix 
Mach number . 



Vll 


P, Q, R Jameson’s upwinding terms 

p pressure 

Q compensation terms 

q magnitude of physical velocity 


Re 

Rf 

Rt 

S 

s 

S to 

u 

u 

u ,v ,w 

u, V, w 

V 


y design 


Vmean 


yt 


radius, radial distance; coefficient of determination 
Reynolds number 
radius of fuselage 
radius of wing tip 

coordinates of the wing’s surface in the auxiliary plane 

wing surface function in the computational domain 

arc length along approximated wake location 

velocity at the edge of the boundary layer 

velocity vector in Cartesian coordinates 

velocity components in Cartesian coordinates 

contravariant velocity components 

velocity vector in computational space 

Cartesian coordinates 

ordinate of the design section 

mean ordinate of the target section 

ordinate of airfoil at trailing edge 


a 


x -r i& 

angle of attack 


* 0 / 



Vlll 


p 


A 
A / 
A t 
5 
Sr 


S' 


C 

7 

r 

V 

4 > 

$ 

p 

sinh 

r 


n 


T2 


5 


angle between the wall shear line and the external streamline of the bound 
ary layer 

transformed boundary layer displacement thickness 
magnitude of change in the airfoil surface in physical coordinates 
user specified trailing edge thickness in units of chord fraction 
central-difference operator defined in Eq. (2-20) 
relofting correction 

boundary layer displacement thickness 
degree of extrapolation coefficient 
ratio of specific heats 
circulation 

flow curvature at the approximate wake location 

averaging operator defined in Eq. (2-20) 

vector differential operator 

reduced velocity potential function 

velocity potential function 

density 

smoothing operator, standard deviation 

hyperbolic sine 

airfoil section thickness 

original airfoil thickness 

airfoil thickness at different Mach number 

momentum thickness 



e 


degree of smoothing coefficient 
coordinates in auxiliary plane 
transformed coordinates 


£ > V 

(, v , C 

Subscripts 


avg 

average quantity 

idle 

forward direct-inverse interface 

idte 

aft inverse-direct interface 

II 

index increment 

h j, k 

grid locations in the £, tj , ( directions 

ky 

value at the wing's surface 

l 

lower surface 

le 

leading edge 

0 

stagnation conditions 

s 

singular line location 

T 

iteration time level 

te 

trailing edge 

u 

upper surface 

w 

wake 

V, z 

components in the x, y , z directions 

oo 

freestream conditions 

*?, C 

components in the ( directions 


Superscripts 

n iteration time level 

o degrees 


^ 03 



X 


TABLE OF CONTENTS 

CHAPTER Pa 8 e 

I INTRODUCTION 1 

II DESCRIPTION OF TAWFIVE 10 

II. 1 FLO-30 10 

11. 2 Grid Geometry 22 

11. 3 Boundary Conditions 36 

11. 4 Boundary Layer Scheme 40 

11. 5 Comparison to Experiment 47 

III INVERSE DESIGN METHOD 49 

III. 1 Inverse Boundary Condition 49 

111. 2 Integration of the Flow Tangency Boundary Condition 51 

111. 3 Relofting 56 

IV REMEDYING SPANWISE INSTABILITIES 59 

IV. 1 Spanwise Oscillations 59 

IV. 2 Success 64 

V RESULTS AND DISCUSSION 92 

V.l Boundary Layer and Wake Effects 92 

V.2 Spanwise Grid Skewness 108 

V.3 Wing Planform Effects Ill 

V.4 Initial Profile Effects 137 

V.5 Pressure Distribution Compatibility 149 

V.6 Grid Refinement Effects 158 

V.7 Fixed Trailing Edge Design 161 

VI CONCLUSIONS AND RECOMMENDATIONS 166 

REFERENCES 168 

APPENDIX 

A DERIVATION OF THE FULL POTENTIAL EQUATION IN CURVI- 
LINEAR COORDINATES 175 

B DERIVATION OF THE Cp EQUATION 179 

XOH 


appendix Page 

VITA 180 



LIST OF TABLES 


Table Pa ge 

1 Results from the analysis of the wings designed with different viscous 

interaction assumptions at a M = .8 and a Re = 24 x 10 6 99 

2 Comparison of the total and wing lift coefficient obtained from a fully 

viscous analysis of the wings designed using different viscous interaction 
assumptions at a M = .8 and a Re = 24 x 10 6 Ill 



xm 


LIST OF FIGURES 


Figure 


Page 


1 Typical examples of overlapping and non-adjacent design regions . . 

2 Staggered box finite-volume cell 

3 Staggered box finite-volume cell with compensation terms defined .... 

4 Three dimensional staggered box finite- volume cell 

5 Conformal grid topography 

6 The computational domain 25 

7 Section surface and wake representation at a constant f station m the 

normalized plane 

8 Section surface and wake representation at a constant r station in the 

auxiliary plane 

9 Constant f surface in the auxiliary plane 

10 Stretching function for the ( (I) direction 

11 Stretching function for the rj (J) direction 

12 Stretching function for the ( (K) direction 

13 Comparison between stretched and unstretched 0 

14 Comparison between the ghost-point potentials defined with flow tan- 

gency (dashed lines) and extrapolation (solid lines) 

15 Comparison of experimental and analytical pressures for RAE Wing- A . 

16 The effect of relofting on the design in the initial stages of convergence . 


6 

14 

16 

18 

23 

26 

28 

29 

30 

32 

33 

34 

35 

38 

48 

57 


17 Alternating thick-thin sections for a divergent medium grid case 

18 The time history of the residual and the terms composing it at the upper 

trailing edge for a typical divergent solution. (The design region extends 
from 30-70% semispan) 


^07 



XIV 


Figure Page 

19 The potential field on the wing’s surface for a diverging design solution . 67 

20 Comparison of airfoil sections designed using the four spanwise oscillation 

remedies on a medium grid from 30% to 70% semispan 72 

21 Breakdown of the Residual at the trailing edge on the upper surface for 

the four spanwise oscillation remedies after 10 global iterations 73 

22 Comparison of airfoil sections designed using the four spanwise oscillation 

remedies on a medium grid from 0% to 100% semispan 76 

23 Comparison of airfoil sections designed using the four spanwise oscillation 

remedies on a fine grid from 0% to 100% semispan 79 

24 The effect of the four spanwise oscillation remedies on the coefficient of 

determination for the medium-grid, partial-wing design 83 

25 The effect of the four spanwise oscillation remedies on the average percent 

error between the target and the design sections for the medium-grid, 
partial-wing design 84 

26 The effect of the four spanwise oscillation remedies on the coefficient of 

determination for the medium-grid, full- wing design 85 

27 The effect of the four spanwise oscillation remedies on the average percent 

error between the target and the design sections for the medium-grid, full- 
wing design 86 

28 The effect of the four spanwise oscillation remedies on the coefficient of 

determination for the fine grid, full wing design 87 

29 The effect of the four spanwise oscillation remedies on the average percent 

error between the target and the design sections for the fine-grid, full- wing 
design 88 

30 Grid System at the wing-body’s surface with a radial boundary stretched 

twice as far as the original grid system 90 

31 Comparison of the unrelofted sectional shapes designed using different vis- 
cous interaction assumptions with a pressure distribution obtained from 
a fully viscous analysis of Lockheed Wing-A at a Re = 24 X 10 6 ,Af = 

.8, and an a = 2° 95 

ao8 



32 Distribution of normalized boundary layer displacement thicknesses on 

the upper and lower surfaces of Wing- A 

33 Comparison of the relofted sectional shapes designed using different vis- 

cous interaction assumptions with a pressure distribution obtained from 
a fully viscous analysis of Lockheed Wing-A at a Re = 24x10 , M = .8, 
and an a = 

34 Comparison of pressure distributions obtained from a viscous analysis of 

the Lockheed Wing-A using three different viscous interaction assump- 
tions with a M = .8, oc — 2°, Re = 24x10 

35 Comparison of an inviscidlv designed section using viscous pressures with 
the target and the target plus the associated boundary layer displacement 

thicknesses 

36 Comparison of pressure distributions obtained by fully viscous analyses 

of the wings designed using different viscous interaction assumptions for 
a M = .8, a = 2°, Re = 24xl0 6 

37 Comparison of the relofted sectional shapes designed using different vis- 
cous interaction assumptions with a pressure distribution obtained from 
a fully viscous analysis of Lockheed Wing-A at a Re = 24 x 10 , M - .8, 
and an a = 2°. The design region extended from 30%-70% semispan. . . 

38 Comparison of pressure distributions obtained from a viscous analysis of 

the Lockheed Wing-A using three different viscous interaction assump- 
tions with a M = .8, a = 2°, Re = 24 x 10 

39 Comparison between a fairly nonskewed (a) and skewed grid (b) 

40 Sections designed with a skewed grid using pressures obtained from an 

analysis on a nonskewed grid 

41 Comparison of a pressure distribution obtained from an analysis of a 
skewed grid with one obtained from an analysis of a nonskewed grid . . . 

42 Grid generated about Wing-C with an incompatible root section and fuse- 

lage cross section 

43 Surface velocity vectors for Lockheed Wing-A on the upper and lower 

surface 



Figure Page 

44 Surface velocity vectors for Lockheed Wing-B on the upper and lower 

surface 118 

45 Surface velocity vectors for Lockheed Wing-C on the upper and lower 

surface 119 

46 Cross-section of the flow in the x-y plane for Lockheed Wing-C 120 

47 Pressure contour plot for Lockheed Wing-A M = .8, a — 2° 121 

48 Pressure contour plot for Lockheed Wing-B M = .8, a = 2° 125 

49 Pressure contour plot for Lockheed Wing-C M = .796, a = 4° 129 

50 Comparison of the designed sections with the targets and the initial sec- 
tions for a fine grid case using Lockheed Wing-A 133 

51 Comparison of the designed sections with the targets and the initial sec- 

tions for a fine grid case using Lockheed Wing-B and a design region 
beginning at 2.5% aft of the leading edge 135 

52 Comparison of the designed section with the target for an unrelofted fine 

grid case using Lockheed Wing-C 138 

53 Comparison of sections designed using two different initial sections .... 141 

54 Comparison of the designed sections with the targets and the initial sec- 

tions for a fine grid case using Lockheed Wing-B ?.nd a design region 
beginning at 5.0% aft of the leading edge 143 

55 Comparison of the pressure distributions obtained from an analysis of the 

Wing-B design, which had a design region that began 2.5% aft of the 
leading edge, with the target pressure distributions 145 

56 Comparison of the pressure distributions obtained from an analysis of the 

Wing-B design, which had a design region that began 5.0% aft of the 
leading edge, with the target pressure distributions 147 

57 Comparison of the designed sections with the targets for a case whose 

design region began 1% aft of the leading edge 150 


Figure 


Page 


58 Comparison of the section designed at a M = .2, using input pressure 
distributions obtained from an analysis at a M = .1, with the original 
NACA 0012 sections 1^2 


59 Comparison of the sections designed at a M = .85, using input pressure 

distributions obtained from an analysis of Lockheed "W ing- A with N AC A 
0012 airfoils at a M = .8, with the original NACA 0012 sections and initial 
NACA 0006 sections ^ 4 

60 Comparison of the pressure distributions obtained from an analysis at a 
M = .85 of the Lockheed Wing-A which was designed at a M = .85 using 
input pressure distributions obtained from an analysis at a M = .8, with 

the input pressure distributions 1^6 

61 Comparison of the sections obtained from a medium grid design at 3VI .8 

using input pressure distributions obtained from a fine grid analysis of 
Lockheed Wing-A with the target sections 159 

62 Comparison of fine grid pressure distributions with medium grid pressure 

distributions obtained from an analysis of Lockheed Wing-A at various 
angles of attack and Mach numbers 160 

63 Comparison of the sectional shapes designed with a fixed trailing edge 

region with a NACA 0012 and the original Lockheed Wing-A section . . 162 

64 Comparison of the pressure distributions obtained by an inviscid analysis 
of the sectional shapes designed with a fixed trailing edge region with those 
for a NACA 0012 and the original Lockheed Wing-A section. (M = .8, 

a = 2°) 164 



1 


CHAPTER I 
INTRODUCTION 

With the advent of efficient numerical schemes that accurately model the ir- 
rotational transonic flow about complex configurations such as wing-bodies and the 
appearance of computers with memory capacities and computational speeds neces- 
sary to execute these schemes in a reasonable amount of time, the efficient design of 
wings for transonic flight is quickly becoming a reality. Although transonic potential 
schemes combined with integral boundary layer solvers may not model the real flow- 
field as accurately as Euler or Navier Stokes schemes , 1-3 their use can significantly 
reduce the costs and time expenditures associated with transonic wing design. 

There are basically two general types of inverse design methods: inverse solvers 
and predictor/corrector (P/C) methods. In the P/C type methods, an analysis code 
is used to calculate the flowfield for an arbitrary initial geometry; and then, this geom- 
etry is systematically modified by considering the differences between the calculated 
and target pressures. The changes to the airfoil sections can be obtained through 
optimization type procedures; or, as shown by Campbell , 4 the appropriate geometry 
changes can be systematically determined by using a design algorithm which relates 
pressure changes to changes in airfoil curvature. 

An example of an inverse solver is the direct-inverse transonic wing analysis- 
design method, which has been under development at Texas A&M University . 5-15 In 
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this method, the wing geometry is determined by specifying pressure distributions 
over part of the wing and then solving the mixed Neumann and Dirichlet boundary 
value problem associated with the full potential equation for compressible flow via 
finite difference and/or finite-volume techniques. The specified pressure distributions 
can be selected by the experienced designer to have such desirable characteristics 
as weak or nonexistent shock waves, a slowly increasing adverse pressure gradient to 
limit boundary layer separation, a center of pressure location giving a desirable pitch- 
ing moment, or an efficient spanwise loading. The designer may also use wind-tunnel 
tests of successful airfoils as an aid in picking a desirable pressure distribution. The 


direct-inverse technique has been successfully used in stretched and sheared Carte- 
sian coordinate systems 5 " 12 ’ 16 * 17 and most recently by Gaily 13 " 10 in a curvilinear 

coordinate system. 

It would be convenient if only the inviscid flowfield had to be included m the 


design process; but, unfortunately, it has been verified through transonic wind tunnel 
tests at low Reynold's numbers and flight testing at high Reynold’s numbers that vis- 
cous effects are very significant 18 . For example, as the Reynold’s number increases, 
the shock wave location is further aft on the wing. Thus, the shock wave m a viscous 
flowfield (finite Re) is located further upstream than that predicted by an inviscid 
(infinite Re) flowfield calculation. Although the inclusion of the viscous interaction 
significantly weakens the shock strength compared to inviscid results, the accompa- 
nying upstream displacement of the shock wave causes the sum of the differences 
between the upper and lower surface pressure distributions to be smaller than in 
the inviscid case; hence, the wing lift coefficient will be smaller in the viscous case. 


Hi 
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Furthermore, it has been discovered that a wing using an aft-cambered airfoil section 
designed inviscidly for transonic conditions might develop 25-50% less lift in a viscous 
environment . 

In light of the previous discussion its obvious that viscous effects must be taken 
into account through some means. One approach that applies in cases where there 
are no regions of massive separation is referred to as the weak viscous interaction 
technique. Since the weak primary viscous interaction effect is the formation of a 
boundary layer on the wing which effectively makes the airfoil thicker, the external 
streamlines for the wing boundary of the inviscid potential field are shifted outwards 
by a distance called the displacement thickness. This shifting is due to the decrease 
in velocity of the fluid in the boundary layer 19 . Thus, to include the effects of weak 
viscous interaction in an analysis of a wing, one simply needs to determine the po- 
tential solution for the surface, find the displacement thickness using the properties 
associated with the streamline representing the body, add this displacement thickness 
to the original surface, and repeat the process until the displacement thicknesses and 
the potential field converge. 

Weak viscous interaction can be included in the inverse design process in much 
the same way. In the inverse regions, where the pressure boundary condition is 
applied, the new surface which approximately satisfies the boundary condition is 
calculated periodically by an integration of the flow boundary condition. At that time, 
the displacement thickness from the boundary layer calculations can be subtracted 
from this new surface to yield the hard or actual designed airfoil. This process can 
be carried out iteratively until there is an insignificant change in the displacements 
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due to boundary layer interaction and the inverse boundary condition, and in the 
flowfield’s potential solution. 

Fortunately, there is a computer program called TAWFIVE (for Transonic Anal- 
ysis of a Wing And Fuselage and Interacted Viscous Effects) which not only has the 
capability of computing the potential field about a wing and fuselage combination 
but also contains a robust three dimensional integral boundary layer scheme which 
provides the necessary viscous effects in the form of boundary layer displacement 
thickness, wake curvature, and wake thickness. It should be noted that a three 
dimensional boundary layer code is desirable in order to properly predict the in- 
creased decambering of the sections near the tip due' to the cross flow in the bound- 
ary layer 20 . In TAWFIVE, the inviscid numerical scheme is based upon Jameson 
and Caughey’s FLO-30 conservative, finite-volume, full-potential flow method where 
computations are performed on a body-fitted, sheared, parabolic, wind-tunnel type 
coordinate system. The three dimensional boundary layer scheme added by Streett 20 
to the originally-inviscid code computes the first order, weak, self-consistent, viscous 
interactions which include the boundary layer displacement effect on the wing's sur- 
face, the displacement in the wake, and the curvature/pressure jump in the wake. 
The boundary layer on the wing is found using a compressible integral method for 
laminar and turbulent flow with a fixed transition location. The turbulent method 
was based on work by Smith 21 , while the laminar method was developed by Stock 22 . 
Small regions of separation are also modeled. This latter feature is an important ad- 
dition for successful convergence, since small regions of separation often occur in the 
initial stages of computations behind shockwaves, in the cove region of aft-cambered 

ais' 



airfoils and near the trailing edge on the upper surface of the wing, even though 
the} may not exist in the final converged solution**. The parameters in the wake re- 
gion are computed in streamwise strips using a two dimensional entrainment integral 
technique. This method has been deemed valid for transport type wings 20 . 

Gaily 13- * 5 has successfully incorporated the inverse design process into the 
TAWFIVE program. Since the modifications made were compatible with the existing 
computational methods and program structure of TAWFIVE, his work resulted in 
a versatile design code capable of allowing the user to design an entirely new wing 
or even discontinuous, nonadjacent segments of a wing. The latter option may be 
invaluable to engineers who are typically faced with the dilemma of designing around 
regions where the wing geometry may be fixed by constraints other than aerodynamic 
considerations. As seen in Fig. 1 these segments can even be non-adjacent upper or 
lower surfaces with overlapping lower or upper surfaces respectively. 

On the other hand, as a consequence of the inverse method, previous experience 
has revealed that specified pressure distributions may not be imposed in regions less 
than about ten percent behind the leading edge of the wing section. This limitation 
was due to the difficulties associated with enforcing the pressure boundary condition 
near the leading edge of the airfoil where the vertical velocities are large. However 
this feature was not viewed as a real limitation since the leading edge regions of most 
airfoils are similar, the leading edge shapes may be constrained by non-aerodynamic 
factors, and since a leading edge geometry can be selected to produce the desired 
pressure values at the beginning of the design region^. 

(p 



Fig. 1 Typical examples of overlapping and non-adjacent design regions 











Moreover, the imposed pressure distributions may often lead to an impractical 
airfoil that has an excessively blunt trailing edge or one in which the upper or lower 
surfaces cross prior to the trailing edge resulting in a fish tail shape. An excessively 
blunt trailing edge might cause a wing to have an excessive amount of drag due to 
base pressure at the trailing edge, while the fish tail shape would be impossible to 
construct. Since the nose shape or curvature has been shown to control trailing edge 
closure. 1 these undesirable shapes can be eliminated with a procedure which 

systematically modifies the leading edge thickness distribution called relofting. Two 
types of relofting procedures have already been included in the program by Gaily. 
One is a simple linear rotation scheme where the surface being designed is rotated 
about the leading edge a proper amount to achieve the desired trailing edge thickness. 
In the second procedure, the leading edge is proportionally thinned or thickened a 

proper amount so that the relofted leading edges are in the same family of airfoil 
shapes. 

Gaily s original design code has been tested in a variety of ways for a Lockheed 
Wing- A wing-bod}. The self-consistency of the approach was tested by designing 
airfoil sections using certain desired pressure distributions, analyzing the resulting 
designed airfoils, and then comparing the desired pressure distributions input to those 
found through analysis. In all of the inviscid cases considered, the code proved itself 
consistent; the section lift coefficients of the designed and target sections and the 
respective pressure distributions were in strong agreement. The relofting procedures 
and the ability of the code to make large surface changes was verified by transforming 
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a 12% thick airfoil at supercritical conditions to a 6% thick airfoil at subcritical 
conditions in the same NACA family. 

Although the code worked well for the inviscid cases attempted, there were 
some modifications and test cases which were required to make this code more valu- 
able. For instance, since Streett found that the wake effects (wake displacement and 
curvature) were relatively important in the calculation of the lift distribution on a 
three dimensional wing, 20 presumably their includsion in the design process would be 
important as well. This was investigated by utilizing the wake options in the code and 
and comparing their effect on the design of a wing. The logic necessary to include 
the viscous effects in the design process originally added by Gaily was tested and 
modified where necessary. 

Recently, a spanwise decoupling in the design regions which led to instabilities 
in the design solution was observed. The supposed source of this instability and the 
various methods used to combat this problem will be discussed later in the report. 

One modification added to the program, which helps smooth out the rippling 
spanwise variations in the wing and give the designer added versatility, is an option 
where the user specifies pressure distributions at the edges of the design region and 
then the changes in the thicknesses of the airfoil sections calculated by the program for 
those stations are interpolated and added to the stations delimited by the edges. This 
approach is different from the original method where the target pressure distributions, 
not the change in thicknesses, were interpolated to the stations in the design region. 

Since the designer is admonished in the TAWFIVE user’s manual 25 that the 
wing is not modeled accurately enough to allow' analysis of very low-aspect ratio 
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wings and that grid problems may be encountered for wings which have high taper 
ratios or sweep angles, three wings of different aspect ratios and sweep angles will be 
used in the inverse design process to approximately delimit the range of geometries 
applicable to the present design code, TAW5D. 

Because of the high computer costs associated with executing this program for 
fine computational grids, results will be shown which will reveal how fine the grid 
needs to be for satisfactory preliminary designs. 

In summary, this thesis presents developments in the inverse design method. 
It includes a brief description of the analysis and design methods and techniques 
used to suppress a spanwise oscillation problem resulting from the interaction of the 
design method with the potential solver. In addition, it presents a series of test cases 
that reveal the lack of dependency of the design on the initial airfoil section, the 
importance of including viscous effects in wing design, and constraints due to aspect 
ratio, wing sweep, spanwise grid skewness. In addition, some questions about the 
necessary refinement of the grid and about any necessary constraints due to Mach- 
number-input-pressure-distribution compatibility will be answered. 
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CHAPTER II 

DESCRIPTION OF TAWFIVE 

As was stated in the introduction, the inverse- wing-design program, TAW5D, 
which was originally modified by Gaily, 13-15 uses as its core the computer program 
TAWFIVE, which can be broken into three major sections: the inviscid, transonic, 
potential flow solver; the cylindrical/wind-tunnel type grid generation scheme; and 
the three dimensional, laminar and turbulent, integral boundary layer code included 
by Streett 20 which is based on the works of Smith 21 , Stock* 2 and Green. * 6- * 8 Since 
the theory behind the code is spread across numerous references, an attempt will be 
made to summarize its formulation in a succinct fashion for the reader s convenience. 

II.1 FLO-30 

The transonic potential flow solver, FLO-30, 29-35 by Jameson and Caughey, 
is a finite volume method which solves the full potential equation in divergence form 

(P u )z + (P V )y + (P w )z = 0 (2-1) 

transformed from Cartesian to curvilinear coordinates : 

(phU)^{phV) n ^{phW\ = 0 (2-2) 

The derivation of the transformation of Eq. (2-2) is presented in Appendix A. 

An expression for the local density, p, and the local speed of sound, a, nondi- 
mensionalized by the appropriate freestream quantities can be found by beginning 
with the energy equation 



l 


(2-3) 



where q 2 = (u 2 + v 2 + w 2 ) 

Then assuming the existence of a perfect gas such that 


h = c p T = 


7-1 


(2-4) 


the energy equation becomes 


7 ~ 1 


7 7 

-<?i + 


7-1 


<?2 4 * a ? 


(2-5) 


Next, assuming freestream and stagnation conditions such that 

<?i = 9oo a i = doc 

<?2 = 0 a7 = a 0 

and upon normalizing all the primitive variables by the appropriate freestream 
tities 


( 2 - 6 ) 


quan- 


f-*r 

pocq oc 
a 


p = 


Poo 


a = 


9oc 


T = 


( 2 - 


The bars on the nondimensionalized quantities will hereafter be omitted for 


conve- 


nience. 


Eq. (2-5) becomes 


; _ 7-1 , 1 


ML 


( 2 - 8 ) 


The local speed of sound is obtained using Eqs. (2-5) and (2-8), yielding 

7-1" 


a = a; - 


Qoo 


(2-9) 


L'sing the isentropic relation 
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and realizing that 


Poo = 


7 Mle 


the isentropic relation becomes 


_ P_ 

V 7 Ml, 

Then making use of the speed of sound relation 

2 7 V 

a — — 

P 

a relation for density is found 


P = (aMoo)^ 1 


which for air can be simplified to 


P = I “ 

^ DO 


7 

7-1 


a 

&OQ 


( 2 - 11 ) 


(2-12) 


(2-13) 


(2-14) 


(2 - 15) 

This expansion is the actual form used in FLO-30, but the more familiar formula for 
density is shown in Eq. (2-16) and can be easily determined by substituting the speed 
of sound relation of Eq. (2-9) into Eq. (2-14). 


1 



(2-16) 


The nonconservative form of Eq. (2*1) shown in Eq. (2-17) can be determined 
by expanding the derivatives of Eq. (2-1); substituting in the appropriate derivatives 
of the density using the expression in Eq. (2-3); multiplying by and then imple- 
menting the equation of state for a perfect gas, the definition of the speed of sound; 

and finally defining the velocities in terms of a velocity potential. 6. 

(a 2 — u 2 ) 6 XZ -f ( a 2 — v 2 ) 6 yv -r (a 2 — w 2 ) 6 ZZ 

(2-17) 

— 2 uVO X y — 2VVJ0y Z — 2 XLWQzz = 0 
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Both of these forms are valid for isentropic, irrotational flows of Mach numbers 
ranging from zero to transonic 25 ; but, by using the conservative form of the potential 
equation, a finite difference scheme will result 36 which conserves mass, especially in 
areas containing large gradients such as with the flow through a shock. Although, 
nonconservative schemes have been successively implemented due in part to the fact 
that the effective mass production at the base of the shock wave fortuitously models 
the shock/boundary layer interactions, the best approach may be to use a conservative 
scheme with viscous corrections added by a separate boundary layer model 3 '. This 
approach is the method utilized by TAW FIVE to include viscous effects. 

FLO-30 uses a finite- volume type scheme which makes use of a staggered box 
approach. Its formulation is directly analagous to the control volume approach used 
to derive the original PDE in Eq. (2-1), except in the finite-volume scheme, the 
discrete nature of the finite difference model is considered from the onset bv usin^ 

* O 

a finite control volume in the neighborhood of a grid point in the finite-difference 
mesh 36 . This method is best illustrated by using it to discretize the following two- 
dimensional. incompressible version of Eq. (2-1) written in Cartesian coordinates 

(2 - 18) 

With the aid of the two-dimensional box shown in Fig. 2, it can be seen that the 
staggered box scheme derives its name from the way in which the primary and sec- 
ondary boxes interlock. The values of the potentials at the four grid points which 
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V 

U+V2 



make up the corners of each primary box are used to calculate the velocities, u . v , in 


the following manner: 


u — o~ =[iy6 z 6 

( 2 - 19 ) 

U — 2T ^ y® 

where fi and S are averaging and differentiating operators respectively and are defined 


by Jameson as 


P-zf 2 (^i+j ,j ' 


(2-20) 

where it is assumed that Ax = 1 . Therefore, the velocity, u. for instance, at the 


primary box center located at (* + + 7 ) is found by 




(©t+l,j ~ + (dj*l,j-M ~ ^IJ-1 ) 


(2-21) 
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The flux at the midpoint of each secondary box is determined by averaging the ve- 
locities u and v at the corners of that box in the y and i direction respectively; and 
the net flux into the secondary box at (i t j) is obtained, giving the discretized version 
of Eq. (2-18) 

Hy&x (u) + fi x S y (v) = 0 (2-22) 


where for example 


{HyS x u) itj = 


V- 11— U l ■ 1 -f V • . 1 - . l — U- l 1 ) 


(2-23) 


The previous discussion implicitly assumes that the velocity varies in a linear fashion 
between the primary cell centers so that the flux into the top of the secondary cell 


face would be, for instance: 
v(x.y)dx 


l 


, 1 , . 2 1 , . 1 V- . 1 ■ . 1 — V; 1 - . 1 




1 , ■ 1 


Ax 


X *r 1 ■ . 1 ax 


(2 - 24) 


t’ ■ 1 _• , 1 “T v • , 1 . , 1 


Jameson and Caughey found that this lumping of the fluxes at the primary cell centers 
reduced to a rotated Laplacian type difference scheme and hence to an uncoupling 
of the solution between adjacent grid points. Therefore, compensation terms were 
added which basically extrapolate the fluxes from the corners of the secondary cell 
to a distance, 6, towards the midpoint of each secondary cell face. Considering Fig. 3 
and using an e ~ .25, the flux, u, at the corresponding grid location (i -f \,j -r |) is 


U- 1 




, 1 = 1 J . 1 “ .25 




du 


dy J i- 1,-i 

1 ' 2 ' 2 


(2-25) 
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Fig. 3 Staggered box finite-volume cell with compensation terms defined 


where 


du 


= <•»(?)(+}. i+s 


(2-26) 


dy / i • 1 ■ 1 

When all the fluxes are extrapolated in this manner, the fluxes at the secondary cell 
centers become 

- £ (Px y )^i^i + U i+ ij_i 4- e(6 X y) i+ i^_i 

2 

" c (©xy) t -l 


u, ' + 5-i 


U- 1 =■ 


r '*^ =• 


V *0-J 


2 


(2-27) 


'XZl 



17 


When the net flux into the secondary box is accounted for, Eq. (2-22) becomes 

p y 8 x u -I- p x 8 y v 

(2 - 28) 

- e [i^x ») t - + i j + » - (^*y),- + » t ;_» ~ ( < £*y)t-}j+$ + (^*y)i- \j-\) ~ 0 
which is equivalent to 

Hyy8xx 4* d“ Pxx8yy<P.~ e ^**yy^ = 0 (2 — 29) 

(Typically, e is .25) 

Notice that the compensation terms lead to a fourth derivative of the potential; this 
higher order derivative will become important later in the discussion of a spanwise 
oscillation problem that occured in the design process. 

The previous concepts can be extended to three dimensional compressible flow 
in curvilinear coordinates by considering eight primary boxes as shown in Fig. 4. The 
three-dimensional potential equation 

{phu\ + (phV) v + {ph\V) c = 0 (2 - 30) 

is again descretized in the same way as in the two-dimensional case to give 

(P hU ) + Hth (■ P hV ) + Huh (■ P hW ) = 0 ( 2 " 31 ) 

The same averaging scheme is used in this case except that the derivatives now have to 
be averaged in two of the coordinate directions instead of one. For example, ( ph\\ 
becomes: 


Me phW)u, k = 




4 


4 



(2 - 32) 




Since relating the potential, <p, to the contravariant velocities, U, V', and, W may 
be somewhat unclear to the uninitiated, it is explained here for convenience. First, 
considering the the full potential function, <$, defined as 

$ = <p + x cos(q) + y sin(a) (2 — 33) 


the standard chain rule can be applied to it to give it, v and,u> as follows: 


<9$ do di do drj do d( 
di dx ' dr] dx d( ' dx 
dodi _ do dr] _ dodi 
di dy ' dq dy ’ d( dy 
do di do dr] do di 
di dz ^ dq dz ' di dz 


U dx 

dy 

_d$_ _ 


-f cos q 




sin a 


dz 


(2 - 34) 


Defining 


(U v= G 

[*7] — £y Vy Cy 
\iz V: G 


(2-35) 


and realizing that 


IA = 


H J 


i-i 


(2 - 36) 


where H is the transformation matrix defined bv 


H = 


x i 

X V 

x c 

Vi 

Vn 

Vi 

Zc 
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ith 


h =|H| 


(2-37) 


the physical velocities, u,v,w normalized by can be related to the gradient of the 
reduced potential function, 6 , by 





-1 f 

\°C/ 



(2 - 38) 


Note that since the grid point coordinate locations in the physical space, (x, y, z), 
are generated as functions of (£, 77 , £), it is convenient to use H instead of J explicitly. 
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The contravariant velocities U,V,W. whose directions lie along the correspond- 
ing grid lines are related to the physical velocities by : 


( U v)=H-'(l 

\wj \v>. 

and the derivatives o{ the potentials and the metrics are defined as: 


(2 - 39) 


<Pv = 


Vi = (») 


(2-40) 


©< = Huh (^) ~ (-) 

The density, ,, and Jacobian, h. are evaluated at the centers of the at their respective 
primary cell centers. Again, by lumping the fluxes at the corners of the secondary 
cell's corners, the solution is decoupled on odd and even grid points leading to two 
independent solutions. This problem is remedied with compensation terms which 
again move the evaluation location of the fluxes to a point somewhere in between the 
corner and the midpoint of the secondary cell face. When this procedure is performed 

for all the cell faces, the potential equation takes on the form of 

fi v< 6 ( ( phU ) + ( P hV ) + hA ( P hW ) 

- e UchvQtv + nhc + Hh&t " 2 6(n(Q(v< ) = ° 

where the Q’s are the compensation terms defined by Jameson as : 

Q( V = ( A ( + A n ) nxkv* 

Qn< = i A v + A () Hh& 

Qtt = (^C + A;) » V S«<? 

Qin C = ( A ( + Ar > + 


(2-41) 


(2-42) 
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Here, A^,A V ,A^ are the influence coefficients which compensate for the dependence 
of p on <p{, 4> v , and <p (. These terms end up being the coefficients of of <f> <f) vv and 

(( in the expanded form of Eq. (2-30) 35 . 

Since the formation of entropy through a shock wave has been neglected through 
the use of the potential function, artificial viscosity must be added to eliminate the 
physically unrealistic solutions. In general, if central differences are used throughout 
the flow field, it is possible for the solution to predict discontinuous expansion shocks 
followed by compression shocks. This situation is a case where entropy decreases 
which is a physical impossibility, and is remedied by adding Jameson’s 30,31 P , Q, and 
R terms which provide the necessary artificial viscosity by producing an upwind bias 
in the supersonic zones. The form of these terms can be found by formulating the 
potential equation in streamline coordinates which reveals the true zone of dependence 
in the supersonic zones. Then in these supersonic zones, the second derivatives of 
the potential, o, included in the streamwise term are formulated with upstream or 
backward differences while the second derivatives included in the crossflow term are 
differenced centrally*®. As shown in the final form of the following finite volume 
equation, the terms are formulated in such a way as to maintain the conservative 
form of the potential equation. 

p v & (phU + P) + p a *n (P hV + <?) + HMP hW + R ) 

( r n r n j c n 6 tvCQ(vC\_n ( 2 ~ 43 ) 

This numerical equation is then embedded into an artificial time dependent 
equation 

A (phU + P)+A(phV + <?)+!: {phW + R) 

< 9 £ or) d(, 

-r compensation terms = *r d* d* 


(2-44) 



and solved via a successive line overrelaxation (SLOR) scheme which sweeps m the 
£ direction along constant < surfaces starting at the root of the wing and implicitly 
solves for the potentials in the 77 direction. Equation (2-43) is a direct statement of 
the conservation of mass and should approach zero as the solution converges. 

After obtaining a solution on a coarse grid, grid halving is used so that the finer 
grid has a better initial approximate solution, thus speeding up the convergence of 
the solution. 

II. 2 Grid Geometry 

The computational grid used by the potential solver, FLO-30, is a body-fitted, 
curvilinear mesh which can be wrapped around a generalized wing-fuselage combina- 
tion that is symmetric about the x-z plane. A body-fitted grid system is desirable in a 
full-potential scheme when the boundary conditions are applied at the actual surface 
of the airfoil. With a body-fitted grid, no interpolation is required and the boundary 
conditions are easily and accurately applied. Because of the shape the grid system 
resembles, it is called a wind-tunnel type grid. An example of this grid is portrayed 
in Fig. 5. The grid shown is the coarsest mesh and has 40 x 6 x 8 points in the £, 77 , 
and C directions respectively. With this grid, the wing becomes a constant 77 surface, 
and each cylindrical looking shell is a constant (,* surface. Constant £ lines can 
seen running spanwise on the wing at constant chord fractions from the leading edge. 
Notice also that due to the conformal transformation used 32 ’ 34 constant £ lines are 
packed close to the leading edge of the wing. This clustering is an attractive feature 
when designing airfoil sections using the direct-inverse approach. Moreover, constant 





25 


£ lines are spaced evenly on the wing and, on the finest mesh, give the designer up 
to 21 spanwise stations where the pressure distributions can be specified. As can be 
also seen from the figure, the lines of constant ( and ij are nearly orthogonal on the 
constant ( surface 34 shown at the wing tip of the airplane, while lines of constant ( 
and < on surfaces of constant V , such as the wing, are not orthogonal except, of course, 
for cases where the wing has no sweep or taper. The lines of constant ( leaving the 
surface of the wing are nearly orthogonal to the surface; this fact will be important 

later on in the discussion of the wing-design methodology. 

The computational grid system is created using a senes of analytically-defined 
algebraic, conformal, and shearing transformations to transform the the wing-fuselage 
combination and surrounding fiowfield in the physical space to a box in the compu- 
tational space shown in Fig. 6. Following Caughey 34 , the polar coordinates r and 6 
are defined in the crossflow planes as 

r = {y~ + -*) * (2 — 4 5) 

6 = tan -1 (2-46) 

The fuselage surface, which is symmetric about the x-v plane, is defined by r — 

Rf(x,8). All points in the fiowfield are then referenced to the surface of the fuselage 

at the same x and 8 location and normalized by the distance between radius, R u of 

the cylindrical surface passing through the wing tip and the radius of the fuselage, 

Rj at the given x and 6 location : 

- _ [ r ~ (2 - 47) 

r ~ [Rt- Rf(x,8)} 
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This normalization causes the lines of constant f, or equivalently k } on the surface 
of the wing to be curved in the x-z plane so they will not coincide exactly with the 
chord line of the airfoil section. This procedure also maps the fuselage to a slit in the 
computational domain. This type of normalization allows for high, low, and mid-wing 
configurations. 

The function Rj(x,6) is found through a Fourier decomposition of the user- 
defined fuselage cross sections such that 

rr, 

R/{xi,0) = Veycosj(* + £) (2 - 48) 

;=i 

The coefficients, a t y\ which are assumed to be continuous functions of x, are spline 
fitted in the x direction for each j. The required radius of the fuselage can be found 
for any point on the wing, or in the flowfield. by interpolating these coefficients to 
the desired x. 

A singular point is located at the focus of a parabola which is fit to the leading 
edge of each wing section with a least squares curve fit. The wing sweep, taper and 
dihedral are accounted for by referencing the coordinates in each surface of constant 
f to the location of the singular line, which is the locus of points comprising the 
singular points, x 5 (f), 8 S (t ) at the leading edge of the wing. 

X = — — + log (2) (2-49) 

c(r) 



(2 - 50) 


} 3 * 
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Fig. 7 Section surface and wake representation at a constant f stauon rn 
normalized plane 


This normalization effectively maps the wing's planform to a rectangle in the compu 
rational space. The , coordinates of the wing corresponding to the given * and X are 
found by linearly interpolating the coordinates of the airfoil sections at input stations 
defining the wing in the spanwise direction. Then at the intersection of a surface of 
constant f with the wing’s surface shown in Fig. 7, the wing section and the wake is 
transformed into a bump in the conformally mapped plane, as shown in Fig. 8, with 
the inverse of the conformal transformation 

x + iB = log [1 - cosh (£ r + tV)] 

a 39 


(2 - 51) 
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Fig. 8 Section surface and wake representation at a constant t station in the 
auxiliary plane 

Fig. 9 reveals an entire constant f surface in the auxiliary plane. A function 5 (f\r) 
is defined to be the t/ coordinate corresponding to the wing's surface defined by the 
input geometry at a constant f . 

The Tj 1 coordinate is sheared out with a simple normalization according to 

( = V = y'/S {{',?) , C = f (2 — o2) 

so that the wing surface lies on a coordinate line in a nearly orthogonal coordinate 
system of £ = const . 

Next, the spacing of the coordinate points in the physical domain is controlled 
by introducing a Cartesian grid into the computational domain where 

— film < £ < Jlim.) 0 < Tj < 1, 0 < C < C Itm 

110 



(2-53) 





Since the derivatives of the spatial coordinates needed for the transformation metrics 
are evaluated numerically, stretching to infinity is impossible; thus the computational 
domain is truncated a finite distance away from the airplane. The outer limits of £ 
and < are chosen such that the grid stretches out far enough from the wing-fuselage so 
that freestream boundary conditions can be safely applied. These constants are not 
user specified, but rather are hard coded in Subroutine COOR of TAWFIVE, such 
that the distance of the outer boundary from the fuselage is about 3 wing spans. This 
distance is probably more than sufficient for most applications; but if a low aspect 
ratio wing is used, which has a large powerful potential vortex at the wing tip and 
significant amounts of spanwise flow, the aerodynamicist may want to increase the 

outer boundary distance. 

The and, C functions for a coarse grid (40x6x8) are shown in Figs. 10-12. 
Notice that distribution of £ between grid points 8 and 24, which corresponds to the 
upper and lower trailing edges respectively, in this domain varies linearly and evenly 
on the wing and then varies quite quickly into the wake ending at a downstream 
location where the flowfield is assumed to be nonchanging. The ^ stretching function 
has the same form, but of course the outer limit at K = 12 determines the outer, 
radial boundary where the freestream conditions are imposed, which in this case, 
as mentioned earlier, will be about 3 wing spans. The t? stretching function varies 
in a parabolic fashion from the wing’s surface at J = 14. Although this stretching 
does seem to pack grid points close to the surface of the wing, since tj is basicalh an 
angular ordinate, the grid spacing above the wing becomes greater as one proceeds 
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Fig. 10 Stretching function for the £ (I) direction 


towards the tip. This increase means that the resolutfon at the tip region is much less 
than that at the root, but this is countered later with a radial correction so that the 
grid spacing immediately above the wing is essentially constant for every spanutse 

station. 

Once the function has been linearly interpolated to the new { coordi- 

nates, the physical coordinates of the grid system can be found through the reverse 
procedure. First, £', V, and f are found using Eq. (2-52). Then Eq. (2-ol) is 
used to extract x and l But before this operation is performed, X and S have to be 



Fig. 11 Stretching function for the 77 (J) direction 


separated in Eq. (2-51), First, both sides are exponentiated and 
hyperbolic cosine is used so that Eq. (2-51) becomes 

A"- i-i(^ + 

Using Euler’s identity. 

e 1 ' = cos(r) -f i sin(r) 


rearranging, and separating imaginary and real parts, gives 


t~ cos 8 = 1 — — cos T\ 



t~ sin 8 = 


1 ■ , 
— sin 77 



the definition of the 


(2-54) 


(2 — 55) 


(2-56) 


(2-57) 
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Fig. 12 Stretching function for the ( (K) direction 


Dividing these two equations by each other and solving for 0 explicit!} }iekL 


8 = tan 


-i 


— sin T)' sinh 


1 — cos 77 ' cosh 


(2 - 58) 


Next x is found explicitly by first using a trigonometric identity and Eq. (2-o8) to 


generate 


6 = 


sin 77 ' sinh £' 


(2 - 59) 


y (1 — cos rf cosh £ f )“ -r sin if sinh £* 
Substituting this into Eq. (2-57) and performing some algebra gives 


= In (cosh - cost?') 


(2 - 60) 


So given C and 77 ' from the previous steps, the normalized coordinates x and 8 are 
obtained for all the grid points in the domain. At this time, two more special stretch- 
ing functions are introduced. One function is used to further stretch x downstream 
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K 

Fig. 13 Comparison between stretched and unstretched 6 


of the wing and another scales 6 such that nearly constant grid spacing is achieved 
immediatel}* above the wing from the root to the outer boundary. The effects of the 
stretching functions can be seen in Fig. 13. 

Notice that this conformal transformation packs grid lines at the leading edge 
of the wing where the gradients are large. This clustering is an attractive feature for 
the inverse design procedure. However, it is paired with the disadvantage that the 
chordwise grid spacing is large at the trailing edge where high resolution is needed to 
accurately satisfy the Kutta condition and to resolve trailing edge pressures accurately 
especially with those generated by aft cambered airfoils. 

Equations (2-49) and (2-50) are inverted to give x and 6 and then Eq. (2-47) is 
inverted to yield r for a given x and 6. This last step requires extensive interpolation 
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to find the radius of the fuselage, R f (x,6), for all of the grid points. Then Eqs. (2- 
45) and (2-46) are used to find the physical coordinates y and c of the grid. Email}, 
coordinates of the points located in ‘ghost ‘ surfaces are obtained through simple linear 
extrapolation of the adjacent grid points along the appropriate or,( S rid line - 

II. 3 Boundary Conditions 

There are a number of boundary conditions which must be applied to the math- 
ematical model of the physical flow about the wing-body. These include flow tangency 
on the wing, fuselage, and the symmetry plane; appropriate far-field boundary con- 
ditions at the finite limits of the computational domain; the Kutta condition at the 
trailing edge of the lifting wing; appropriate treatment of the wake; and the compu- 
tational slit outboard of the wing tip. 

II. 3.1 Flow tangency 

The flow tangency condition is easily implemented due to the curvilinear system. 
The fluxes above the surface need only be reflected to the ghost points beneath it so 
that the net out of plane component of the flux vanishes at the surface. In the case 
of the wing this becomes 

P hU i,ky+\. k = P hU Uky-lk 

phV i ky+ i tfc = -phV iM _ i ifc where: ky = j wing (2-61) 

P hW i,ky+$,k = P hW i,ky~lk 

Similarly for the symmetry plane 

P hU i,\\,k = P hU i,2\,k 

ohV- i , = -ph,V.-i . where: j = 2 on the symmetry plane (2 - 62) 

P ' 1,1 r 

P hW i.l\,k = P hW i,2\,k 



While for the fuselage this becomes 


P hU i,j,2\ = phU i.3' 

phV^^i = phVij 3 1 where: A = 3 on the fuselage (2-63) 

The previously discussed compensation terms and upwinding terms are also similarly 
reflected in an appropriate manner. 

Potentials at the ghost, points located at grid points beneath the surfaces are 
needed for the calculation of surface velocities used in the upwinding terms and the 
surface pressures. These are found for the wing and fuselage by setting the appropriate 
contravariant velocity to zero in 



H J H 




Of 


°C. 


-H‘ 



(2 — 64) 


and using the resulting equation to solve for the unknown potential at the ghost point. 
In the case of the fuselage, this method of defining the ghost points is used solely when 
they are needed in the calculation of the upwinding terms in the residual expression. 
When the pressures are calculated, the ghost points are defined by assuming 


4>a - 0 


(2-65) 


so that the potential at the ghost point is, in effect, linearly extrapolated in the span- 
wise direction. As seen in Fig. 14, these two methods lead to quite different values. 
The first leads to a discontinuous spanwise variation in the potential while the second 
has a much smoother variation. The first approach guarantees that the flow will be 
tangent at the fuselage, while the second does not. However, the pressures calculated 
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at the root are fairly independent of the method used to define the potentials at the 
ghost points in the fuselage. 

The potentials at the ghost points of the symmetry plane are similarly calculated 
by assuming 

d„ = 0 (2-66) 

This process imposes an inflexion point on the pertubation velocity in the 77 direction 
at the symmetry plane since only symmetrical cases are treated. It is uncertain why 
g-q was not set to zero instead to approximate tangency at the plane of symmetry. 
However, this situation is rather academic since these ghost points are \rsed only for 
supersonic regions adjacent to the symmetry plane to compute the small spanwise 
upwinding term. 

11. 3 . 2 Far-field boundary conditions 

Since the reduced potential used in the formulation of the numerical method 
represents a pertubation from the freestream value, the3 r are set explicitly to zero on 
the radial boundary, (max> and the upstream boundary represented by part of the 
minimum 77 surface. 

At the outflow boundary, (£ = £ m j n max)’ streamwise pertubation velocity, 
is set to zero. This latter condition implies that the pressure will return to its 
freestream value, assuming that there is not any crossflow 34 . 

11. 3. 3 Wake treatment 

In the original method of FLO- 30 , the wake is treated as a vortex sheet which 
has a discontinuous jump in the tangential velocity and a continuous normal velocity 

a 
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through the sheet. The rolling up of the sheet is ignored and the vertical convection 
of the sheet is approximated by assuming that the wake lies along the constant rj 
grid line that leaves the trailing edge smoothly and returns to the plane of the wing 
at the outflow boundary. The requirement that the normal velocity be continuous is 
enforced by setting V', = 0 on the wake, which fixes the values of the potentials at the 
ghost points, and the jump in the tangential velocity is satisfied by forcing a constant 
jump in potentials on the the surface of the sheet along a constant ( and t] line. This 
jump in potential is obtained using the circulation determined at the trailing edge of 

the wing. 

11. 3. 4 Outboard computational slit 

Due to the C-grid type formation of the grid, there exists a computational slit 
outboard of the wing tip on the plane of the wing. Since physically the pressure must 
be continuous across this cut, the potentials on the surface and at the corresponding 
ghost points are defined such that the reduced velocities normal and tangential to to 
the surface are continuous across the slit. 

11. 4 Boundary Layer Scheme 
II.4.1 Integral method 

Streett 20 included an integral boundary layer scheme in TAWFIVE to account 
for the necessary viscous effects in the form of the boundary layer displacement thick- 
ness, wake curvature and wake thickness. An integral method was chosen for its 
computational efficiency and its relative robustness. 
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In an integral approach the degree of the partial differential equations is reduced 
by an a priori integration in the direction normal to the surface. 21 This reduction can 
be illustrated by considering the boundary layer equations governing a two dimen- 


sional incompressible flow 


19. 


du dv 
dx ' dy 


(2-67) 

(2 - 68 ) 

If Eq. (2-68) is integrated with respect to y from the wall (y = 0) to a distance h 


dv du dV d~u 

U ~ T* V — l “7 ~ ^ *> 

ox oy dx oy- 


outside the boundary layer, it becomes 


/ / 0V r ° 
/ (u— -r V— - U — )dy = — 

/ 0 dx oy dx p 


(2 - 69) 


where r D is the shearing stress at the wall. 


Using the continuity equation. Eq. (2-67), to obtain the normal velocity component, 


v , as 


1 Jo dx 


)dy 


(2-70) 


and substituting this result into Eq. (2-69), the result is 

f h . dv dv [ v dv , efU t 0 

/ A~ a y - ^T~) dy = ~T 

J y = o ox dy J o dx dx p 

After integrating by parts and reducing, Eq. (2-71) becomes 


(2-71) 


^ [u (U - u)] dy -f ^jr- [ (U — u) dy = — 

Jo 


(2 - 72) 


[ h 

Jo dx * ' " dx J 0 ' ' p 

Now’, taking h — * oo and defining a displacement thickness, 6 j, and a momentum 
thickness, # as 


= / (U-u)dy 

Jy=0 

roc 

9X)~ = / u(U - v)dy 

jy = o 




(2-73) 


and substituting them into Eq. 2-72, it becomes 


In this reduction process, two partial differential equations have been replaced by 
one ordinary differential equation. Since only the integrated quantities. ** and 8, 
are really the only quantities required of the boundary routine to model the weak 
viscous interaction, the fact that the solution to this equation does not provide the 
exact local variation of primitive flow properties across the boundary layer is not of 
consequence. The required functional form of the variation in v across the boundary 
layer is assumed a priori - by a polynomial for instance. 

11. 4. 2 Laminar scheme 

In three-dimensional, compressible, laminar flow the same integration proce- 
dure is implemented using two bounday-layer momentum equations and their corre- 
sponding moment of momentum relations to yield a system of four coupled partial- 
differential equations. 20 In the formulation of these equations, it is assumed that the 
streamwise velocity profile is of the Faulker-Skan (F-S) family of similarity profiles 
and that the cross flow profile is a linear combination of the F-S familj of profiles. 
These incompressible profiles are extended to compressible flow by the scaling of the 
normal coordinate with the Stewartson transformation. 

11.4.3 Turbulent scheme 

The formulation of the turbulent scheme is similar to the laminar, but the 
streamwise velocity is assumed to have a simple power-law profile which is a function 
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of the streamwise shape factor and the transformed boundary layer thickness and 
normal coordinate; and, the cross flow profile has the form of 

F = TT (1 “ —)* (2 - .o) 

U U Za 

where Z is the transformed normal coordinate, A is the transformed boundary layer 
displacement thickness, and 0 is the angle between the external streamline of the 
potential flow and the wall shear direction. In the turbulent scheme, the final three 
governing equations are two momentum integral equations derived from the conti- 
nuity and boundary layer momentum equations and one entrainment equation. The 
latter equation accounts for the addition of mass into the boundary la)er from the 
surrounding flow as the boundary layer grows. 

II.4.4 Lag entrainment 

Originally, in the work by Smith" 1 , the relationship between the entrainment 
coefficient and the shape factor required in the previous scheme was formulated em- 
pirically with a simple algebraic equation. Later Green found a relationship for the 
required quantities through the use of the turbulent kinetic energy equation which 
explicitly represents the balance between production, advection, diffusion and dissi- 
pation of turbulent energy in the boundary layer. He referred to this as the Lag- 
Entrainment method 2 L 

Also, in Green’s method the desired momentum and displacement thickness of 
the wake is determined by simplj' continuing the integration of the three governing 
equations past the trailing edge on either side of the wake . It is assumed that aft 
of the trailing edge that the skin friction coefficient is zero and that the dissipation 
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length scale is twice that on the wing. Once the integration is performed on ejther side 
of the wake, the required integral properties are simply the sum of those calculated 

on both sides. 

II. 4. 5 Solution of the governing equations 

The resulting governing equations are solved through an explicit u pe integra- 
tion scheme in the x (or chordwise direction) along constant span stations. In this 
scheme, the domain of dependence is conservatively assumed to lie between the ex- 
ternal streamline of the potential flow and the shear angle of the boundary layer. To 
account for this dependency, the span wise derivatives found in the governing equa- 
tions are backward differenced if the external streamline and the wall shear line lie 
on the outboard side of the chordline and central differenced if the streamline and 

the shear line lie on opposite sides of the chordline. 

Boundary conditions are required at all inflow boundaries. At the root, a plane 
of symmetry is assumed. Here, the cross flow velocity is set to zero, as are all all 
spanwise derivatives. At the wing tip, all spanwise derivatives are also set equal to 
zero. And finally, an attachment line approach 38 is used to determine the initial 

conditions at the leading edge. 

IL4.6 Wake curvature 

When the flow leaves the wing at the trailing edge, it initially follows a curved 
path and then soon aligns itself with the freestream downstream of the wing. This 
large curvature of the flow near the trailing edge can have a measurable effect on the 
overall lift of the wing. In fact, Streett found that in one instance the sectional lift 
coefficient near the tip of the wing was decreased by about four percent when the 

a. trS' 



curvature of the wake was taken into account. Usually, if only first order effects are 
considered, the pressures at the trailing edge would be equal on the upper and lower 
surface. But, if the wake is considered to have an effective thickness of 6* - 6 due 
to viscous effects and curvature, the pressures on either side of the wake will not be 
equal except at the centerline of the wake. Since the fiowfield about the wing and the 
wake with the displacement thickness added to it is modeled inviscidl}, the trick is 
to calculate a pressure difference across the wake at the trailing edge in the inviscid 
flow which will yield a zero pressure difference at the centerline of the wake in the 
real viscous flow 39 . It has been shown that the appropriate pressure jump across the 
wake with a thickness of S w can be written as a function of the curvature, of the 
centerline of the actual wake, the mean tangential velocity, u u ., and the mean densit}. 
in the wake as 

&P = Ptop Pbottom KpvjU'u-Qv.' 

Given that the pressure difference is small, this can be related to the circulation, F, 
by 

f ~ dT w = - f 8u.Ku.dSw (2-77) 

Jztt J ~tt 

where S w is the arc distance along the wake. The circulation at the trailing edge 
is calculated by the difference in the potentials at the trailing edge in the inviscid 
solution and Eq. (2-77) is numerically integrated from the trailing edge to one grid 
point upstream of the downstream boundary. The circulation at the downstream 
boundary is then matched to the circulation obtained from the integration. 
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Since the wake effects are relatively small, 39 it is only important to know the 
approximate location of the wake centerline. This simplifies the problem since the 
actual wake location would have to be found by tracking the streamline of the inviscid 
solution leaving the trailing edge and then a new grid would have to be created 
about the new wake so that the boundary conditions on the wake could be applied. 
Alternatively, the approximate shape of the wake can be found by assuming thal the 
streamline leaves the trailing edge smoothly at the average of the local trailing edge 
angles and that then the angle between the wake centerline surface and the freestream 
decays logarithmically, similar to that of a point vortex in a uniform freestream at a 
given angle of attack 20 . The circulation, T. of this point vortex located at the quarter 
chord point could be determined by forcing flow tangency at the trailing edge of the 
wing section. The ordinate of the centerline of the wake would then have a form 

similar to 

3 (wake = Vte + tan a(d - |c) - ^Ctan(a)/n^/d (2 - 78) 

where d is the x distance from the quarter chord point of the wing section. 

The curvature of the flow, k, can be determined by calculating the rate of change of 
the flow angle at the approximated wake location. 

II. 4. 7 Wake thickness 

The thickness of the wake is accounted for by simply adding the displacement 
thicknesses obtained from the boundary layer solver to either side of the predefined 
wake location. The ghost points in the wake are then redefined such that strict flow 
tangency is enforced along this new surface. 

1 
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II. 5 Comparison to Experiment 

TAWFIVE was used to analyze RAE Wing-A wing-body at a Mach number 
of .8, an angle of attack of 2 degrees, and a Reynolds number of 2.66 million based 
on the root chord. The pressure obtained from this analysis are compared to some 
experimental data at two convenient stations in Fig. 15. Even though no attempt 
was made to try and match lift coefficients by changing Mach number or angle of 
attack, the comparison between the experimental and predicted pressures is fairh 
good up to the trailing edge. There TAWFIVE predicts slightly higher pressures. 
This characteristic behavior has been attributed to the improper modeling of the 
the strong viscous-interaction region at the trailing edge* 0 but may also be due to 
a combination of the coarseness of the grid at the trailing edge and vrind tunnel 
interference errors. 





Fine Grid 
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CHAPTER III 

INVERSE DESIGN METHOD 


III. 1 Inverse Boundary Condition 

As stated earlier, in the direct-inverse method a pressure boundary condition 
is enforced rather than flow tangency aft of the portions of the wing which are to be 
designed. Following Gaily, 13-15 the input pressure coefficient can be written in terms 
of the Mach number, M, and the freestream speed, q, as 

1 

(3-1) 


C — 

p 7 M- 


1 *> /*),*), 0 \ *> 
where a/ = \tr — v~ — w~) q~ ^ . 

Solving for v in Eq. (3-1) yields 


■ 

7 — 1 2 / q 1 V 

7 -j 


It 0 Moo' ( 1 o ) 

- 1 


2 V q*/J 



u = 


1 - 


(i-TjltT 3 ! 


(l -i- 1 _ i 


( 3 - 


i+(5)-+«r 

This form of the equation seems to have been chosen over the more obvious form of 

(3-3) 


v = 


1 - 




2 

oo 


•v — 1 


1 + 


'rM 2 oc C J> \ ' 


- 1 


— v — w 


since it is less likely that its radicand would be negative. Equating Eq. (3-2) and the 
first row of Eq. (2-38) results in 


+ J\2<?n + ^13 <?c = 




T-l 


(l + 7 _ i 






— cos a (3 — 4) 
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where J itj are the elements of (h T ) \ A potential, can be formulated 

in terms of the pressure coefficient by expanding about the grid point location 
(;_ i J,k), and then using central differences in the ( and ( direction and second 
order backward differences in the normal direction, 77 , yielding 


Jll ~ tf-lj.k 


J 12 [3 - 4 (©fj-i.ir * 

+ i/ 4 (3_5: 

«^13 (oliMl ®i- lj.k-rl — ~ ^ 


=n c ,^ 

where 6 71 are the potentials at the current time level and o 71 ^ 1 are the updated 
potentials. 

Solving for the O to be specified, Eq. (3-5) becomes 

1 


-TI+1 

?iiky,h 


Jn t |^i 2 
— J\2 


JnPi-l,ky,k 

3&x-l,ky,k ~ 4 {&i,ky- 1.* 


/4 


■5* ©a-j/- 2 ,fc + l,fcy- 2 ,fcj 
— ^13 (^”,Jfcy,k+l 4 " <?:-!, fcj(,k-l ~ Qi'ky.k-l ~ ®l-l ,ky ,k -l) / 4 

+ f ( C p.-§.‘)} 


(3-6) 


where F (C p { _i <h ) is the right hand side of Eq. (2-72) and j = ky on the wing 
surface. Also, the 77 grid lines are numbered such that ky - 1 is the location of the 
grid point immediately above the wing’s surface. Pressures are specified at half grid 
point locations in the £ direction to eliminate the chance of the solution decouplin D 


cl ^ / 



on ‘odd’ and ‘even’ grid points. Since the actual sectional shape of the final wing 
is unknown initially, the potentials are specified on the wing's surface at the present 

time level. 


III. 2 Integration of the Flow Tangency Boundary Condition 

Since the grid is boundary conforming, the wing sections in the design region 
must be updated every so often by integrating the flow tangency condition written 
in curvilinear coordinates. After Gaily, the curvilinear form of the equation can be 
found by first considering the flow tangency condition for Cartesian coordinates 

U T VF = 0 with F(x.y,z) = 0 (3 - 0 


where U is the physical velocities and F is the function describing the surface of the 
wing. 

The physical velocities can be related to the contravariant velocities using the 


aforementioned relations, which are repeated here for convenience. 


[Uj = 



= [Hj[V] = HV 


(3-8) 


By using the chain rule in the same manner in which the above expression was derived, 
the gradient, V, of the surface function, F, with respect to the physical coordinates, 
x,y,z can be related to the gradient, V', of the surface function S{(, 77, C) by 


[VF] = 


U 


G 


iy Tjy Cy 
Is Vs C z 



If = (H -1 ) t V'S 


(3-9) 


Substituting these two into the tangency equation gives 


(HV) r (H _1 ) T V'5 = 0 


(3-10) 



Using the identity from linear algebra, 


[A B) t = B t A t 


(3-11) 


Eq. (3-10) becomes 


V r [H _1 H] T V'5 = 0 


(3 - 12) 


which is reduced to the desired form of the flow tangency condition for curvilinear 


coordinates 


V 7 . V ‘S = 0 


(3-13) 


A more convenient form is obtained by expanding this to 


u— + = 0 

d£ dr] d£ 


Since the wing is a surface of constant tj , where 

S{t,V, C) =VU-.C)ky -V = 0 
dS dvitiOky 


= - 1 


di di 

qS_ 

dr] 

55 _ dr,U,Qky 
d( d( 


Eq. (3-14) reduces to 


oUi 


\ _ V W fchf\ 

c.) w 


(3-14) 


(3-15) 


(3 - 16) 


The integration of this equation can be handled in two different ways. If the 
span wise term, is lagged one global iteration, it will always be zero since upon the 



creation of a new grid, all derivatives of T} with respect to the ( or ( direction vanish 
on the wing’s surface; and, Eq. (3-16) reduces to 



(3 - 17) 


The other approach would be to integrate Eq. (3-16) iteratively. If the contravariant 
velocities are frozen at their current values, and the spanwise terms are initially as- 
sumed to be zero, Eq. (3-17) can be integrated to find the approximate inverse changes 
At], These can be used to find approximations to the spatial spanwise derivative, 
which can then be included in Eq. (3-16) to provide a better approximation to the 
flow tangency equation. The process can then be repeated using Eq. (3-16) until the 
spatial derivatives converge. Numerical experiments reveal that the spanwise terms 
are at least two orders of magnitude smaller than the chordwise terms prior to the 
creation of the new grid. Hence, the spanwise terms can normally be neglected. 
Equation (3-17) was integrated using the trapezoidal rule 


II 

Vi,k = — 


II = -1 
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i,k \ / i-JJ,*, 

upper surface 


II — —1 lower surface 
For comparison purposes the fourth order scheme 
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i,k 
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U ' i-IJ,k V U / i—2ll,k 


(3-18) 


t— 3/J,fc / (3-19) 


+ Mi- 


i-II,k 


was also used. With the fourth order scheme the trapezoidal rule was used for the 
first two integration steps. This higher order integration scheme had little effect on 



the final answers, except for coarse grids in regions of high curvatnre such as the cove 


region of a supercnticftl airfoil. 

Since Gaily found that calculating V using strictly finite differences uas not 
accurate enough, he instead, using an approach similar to that in Ref. (60), discovered 
that V was most accurately obtained from the residual expression. First, assume that 

V^J^ = m ipkY) ( 3 - 20 ) 

U * phU PtodpW) 

and then combine the previously defined averaging and differencing operators 


Pn 


(phv ) u> .t - 1 (wVm+WI'U-iJ 13 - 21) 


W = {l* v kkr-l* - (3 22) 


to generate 


S,(phV) iMJt = 2{phV U^-W^'W < 3 23) 

Substituting this result into the residual expressron, Eq. (2-43), and solving for the 


out of plane flux, ph\' , on the wing surface gives 

2 fi w {phVhky* = PvC S ( (pWhky* + 2 ^c {p hV ito-lk) + ( 3 _ 24 ) 

fi iv 6 c (phWi'ky,k) + compensation and upwinding terms 
Since at convergence the flow should also be tangent to the designed surface, 
the tangency condition is enforced in the residual expression, Eq. (2-43), by setting 


(pM'Xky+lk = -(P hV kky- 




(3-25) 



The resulting expression is identical the RHS of Eq. (3-24), and the expression for 
the normal flux becomes 




Residual 

10 


(3 - 26) 


Note that since the residual is not zero in the design region due to the inverse boundary 
condition, this expression reveals that there will be a mass flux of fluid from the 
boundary 23-3 ' during the iterative design process. No attempt was made to account 
for this transient flux, since at convergence it would be zero. 

Upon substitution of Eq. (3-26) into Eq. (3-20) and using the cell averaged flux, 
phU, on the surface the boundary condition becomes 

dy _ V _ Ptv<;{P hV ) _ Residual ( 3 -^ 7 ) 

di U ~~ p (v< {phU) 2 n M (pkU) 

The changes normal to the surface at each spanwise station are obtained by integrating 
from the beginning of the inverse region to the trailing edge using the trapezoidal rule. 

Assuming that the grid line leaving the wing in the 77 direction is normal to the 
wing, these changes, Ary, are then converted from computational to physical units by 
scaling by transformation metrics such that 

{3 ~ 28) 

After subtracting the boundary layer displacement thickness from the inverse changes, 
A/’s, which are linearly interpolated to the user defined input stations, the resulting 
displacements are added to the initial airfoil sections yielding the new wing surface 


for the current time level. 
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III. 3 Relofting 

Many times the trailing edge thickness may be too large if the leading edge 
curvature is too small or may be ‘fish-tailed' if the leading edge curvature is too 
large. These undesirable situations can be remedied by a procedure called reloftmg 
where the designed surface is rotated about the leading edge to meet a specified 

trailing edge ordinate or trailing edge thickness. 6 " 4 

This relofting procedure can be accomplished in two separate ways. 13,14 In the 
first method, assuming both the upper and lower surfaces of the wing are being 
designed, the user specified trailing edge ordinate, 


^ A; 

yt-uvvjcT. IJavg — - ^ 

1 rrti' t r 


(3 - 29) 


is subtracted from the ordinate of the displaced surface, 


V design — V initial vv£tr_ — A ££££! 

leni' er 


to yield a correction of 


&Tt e — 2/i Vdesign 


(3 - 30) 


(3-31) 


where Aj is the user specified trailing edge thickness, is the initial inverse 

change, y tn itial is the trailing edge ordinate of the original airfoil section, and y avg is 
the average of the trailing edge ordinates of the input geometry. 

This correction, 


6 r (x) 


&r tt X 


X ~ XU 

chord 


(3 - 32) 


is proportionally added to the initial inverse displacements -which amounts to a rota- 
tion of the displaced surface about the leading edge to meet the trailing edge ordinate. 


2. k 7 
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Fig. 16 The effect of relofting on the design in the initial stages of convergence 


To illustrate this relofting procedure, the first global iteration of a typical design be- 
fore and after relofting is revealed in Fig. 16. 

If only the trailing edge thickness is specified, allowing the trailing edge ordinate 
the freedom to vary, the correction instead becomes 



A u t Aj 


2 


chord J 


(3-33) 


where Au and Aj are the initial inverse changes on the upper and lower surfaces 
respectively. It should be noted that the inverse displacements are positive when 
they cause an increase in thickness. 

2.&Q 
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The second relofting scheme determines the displacements aft of the direct- 
inverse junction of the design region in the same way, but the leading edge ordinates 
are thinned to meet the displaced surface at the beginning of the design region. This 
insures that the leading edge shapes remain in the same family of airfoils. 


( 3 - 34 ) 


where y id u is the airfoil thickness at the direct-inverse interface in the chordwise 
direction. 

In order afford the designer extra flexibility, one more relofting scheme was 
devised where a portion of the trailing edge region is user specified instead of just 
the trailing edge ordinate. Using the same rational as with the rotation scheme, the 
correction added to the displaced surface to meet, the specified ordinate at the aft 
direct-inverse junction located at Zidtei 1S 


<5 r (z) — ^ 



(3 - 35 ) 





CHAPTER IV 


REMEDYING SPANWISE INSTABILITIES 


IV. 1 Spanwise Oscillations 

In the original work by Gally^’^, the pressure distributions applied at the 
computational grid stations of constant £ lines on the wing in the design region were 
obtained by spanwise, linear interpolation of the pressures input by the user at de- 
sign stations to every grid station delimited. This meant that the inverse boundary 
condition was enforced at every constant (, grid station in the design region, and that 
every sectional shape was determined relatively independent of the others. Unfortu- 
nately, an annoving divergent spanwise oscillation problem sometimes occurred v hen 
designing a wing which required extensive relofting, especially when the initial section 
was thinner than the target. This oscillation led to sections which were too thick or 
too thin at adjacent constant £ grid station, (see Fig. 17). This problem was more 
pronounced when the sweep was increased or the aspect ratio was decreased and was 
usually divergent except for very high aspect-ratio wings (AR=10) with no sweep. 

Early in the research, it was discovered that the problem could be circumvented 
by specifying the C v distribution at at least every other constant ( grid station and 
then linearly interpolating the inverse displacements calculated at those grid stations 
to the other grid stations included in the design region. The regions in the middle 
of the design region were simply analyzed using the original flow-tangency boundary 
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Fig. 17 Alternating thick-thin sections for a divergent medium grid case 
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condition. The resulting sections, interpolated to the geometry input stations, were 
all relofted as usual to satisfy the trailing edge ordinate condition. This procedure 
led to a convergent solution most of the time, except when designing wings with 
significant sweep or with low aspect ratios, such as Lockheed Wing-B and Wing-C. 

It was later discovered that a similar procedure was briefly discussed in Ref. 40 to 
overcome a decoupling of the solution in the chordwise and spanwise direction leading 
to a numerical instability when using an inverse panel-method code. In this case, the 
ordinates of the ‘odd' points along the chord were obtained by quadratic interpolation 
using the ordinates of adjacent ‘even* chordwise points while the ordinates of each 
‘odd' spanwise grid station were generated using linear interpolation between the 
contiguous ‘even 1 spanwise stations. This procedure effectively eliminated half of the 
unknowns. The similarities of the decoupling problem in this scheme and our direct - 
inverse method are quite evident, even though the design schemes are quite different 
in methodology. 

Although this somewhat heuristic cure to the problem seemed to work for the 
most part, the fundamental cause for this problem was not well understood, hence the 
oscillation problem was investigated in much greater depth. Initially, it was thought 
that either the inverse boundary condition or the relofting scheme was solely to blame, 
which led at first to a series of reformulations; while none of these were successful, 
they did create great insight into the problem. 

Since the oscillation problem seemed to stem from the uncoupling of the solution 
in the spanwise direction, the original inverse boundary condition in Eq. (3-5) was 



rewritten as 


, nJ .j 4Jn + 3Ji2 , n +l _ j.n+1 

9uky,k-\ J u i,kv*+l 


— {^n iy.it 

*M3 t 

- J\2 [3^-U-y.i- - 4 + Jfey-l.k) ^ ” 1 } 

-r G'Uy-’.i + ^"-l.fcv-2.*]/ 4 ~ ^ p »-3« fc ) ) 

such that the <?-5 could be obtained implicitly in the spanwise direction. Although 
this would seem to strongly couple the potential field in the spanwise direction, it did 

not deter the solution from oscillating in the slightest regard. 

One form of Eq. (3-4) was tried using one-sided differences for the spanwise 

derivatives, and yet another which specified the C v at (i -r \,ky,l t - 5 ) S nd locations, 

but thev did not cure the problem either. 

The idea of devising a conservative formulation of the inverse boundary con- 
dition using a control volume approach more in keeping with the spirit of the finite 
volume scheme used in FLO-30 or the approach used in Ref. 


( \ 








the details necessary to implement this approach were never pursued. 

Attention was then directed towards the methods used to integrate the flow 
tangency equation and the relofting of the resulting shapes. Since the problem seemed 
to stem from the lack of spanwise information, the spanwise terms m Eq. (3-16) 
were included during the surface update process. The ratio ^ was obtained from 
Eq. (2-39) and the potentials at the present time level. An approximation of the 
spanwise derivatives, §J, was calculated using central spanwise differences of the 
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initial displacements which were calculated using Eq. (3-17). Then Eq. (3-16) was 
solved iteratively until there was no appreciable change in the displacements. In case 
the relofting adversely affected the results, this process was also tried after the inverse 
displacements were changed with relofting. However, the inclusion of these terms had 
very little effect on the displacements calculated since, in both cases, they were at 
least an order of magnitude smaller and did not help the divergence problem in the 
slightest regard. 

Spanwise smoothing of the displacements was also tried. Although this tech- 
nique did provide a smoothly varying distribution of sectional thicknesses, the diver- 
gence was merely slowed. Sometimes the solution would reach a settling point where 
it would not converge further, but the resulting section shapes were not satisfactorily 
accurate. 

In the midst of the search for a cure for the oscillation problem, it was discovered 
that if the potentials obtained from a converged solution of the target section were 
specified on the wing using a different initial geometrj r , the design solution would 
converge without oscillating. This result appeared to condemn the inverse boundary 
condition and redeem the integration and relofting schemes. On the other hand, if 
the inverse boundary condition was applied at every grid station, and displacements 
were calculated only at every other spanwise grid station and were interpolated to the 
stations in between, the solution also converged, which seemed to indicate that the 
inverse boundary condition was not the sole origin of the problem. Thus, it appeared 
that the problem was stemmed from a combination of causes. 

Zl 1 -} 
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IV. 2 Success 

After the many failed attempts of remedying the oscillation problem by refor- 
mulating the inverse boundary condition and the integration and relofting shemes. 
attention was directed towards the residual and the terms composing it. The residual 
is directly affected by the inverse boundary condition; moreover, the residual direct!} 
influences the section shapes through the integration of the flow tangencv boundary 
condition. Consequently, the residual was broken into its major components and 
plotted in the spanwise direction after each surface update of a known divergent case. 
This case happened to be a medium-grid design of Lockheed Wing- A with the initial 
section being a NACA 0006 section over the entire wing and the target being a NACA 
0012 section. The design region extended from 30% to 70% semispan. Sample plots 
for this divergent case are shown at four different time levels in Fig. 18, where the 
total residual also includes the upwinding terms. As can be seen, the compensation 
terms, which include spanwise derivatives of 6 , at first are very small compared to 
the rest of the terms but later tend to dominate and amplify the oscillation. This 
oscillation starts at the direct-inverse interface or, in other words, at the first span- 
wise station from the root in the design region and propagates spanwise as a damped 
oscillation with a period of two grid spacings. 

The oscillation problem seems to be driven by a combination of events which 
build upon each other causing a divergence. It is believed that the initial mismatch in 
the potentials at the direct-inverse interface in the spanwise direction is amplified by 
the compensation terms which include spanwise derivatives of the potential function. 
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The residual is then undershot and overshot on alternating spanwise stations. This 
oscillation is further magnified by relofting, which creates a section that is too thin 
when the slopes defined in Eq. (3-27), which of course are directly proportional to the 
residual, are too large and vice-versa. Since more or less fluid has to be ejected from 
the section that is too thin or thick, respectively, to give a streamline approximately 
corresponding to the correct target section, the potential field shown in Fig. 19 
at each design station is forced further away from the adjacent fields by the inverse 
boundary condition which in turn forces an even further undershoot or overshoot of 
the residual, resulting in a growing span wise oscillation. With the aid of other nu- 
merical experiments, it has been found that it is only necessary to have two adjacent 
design stations to drive this oscillation to divergence. It is of interest that when the 
wavy wing surface resulting from a divergent solution was analyzed with TAW FI\ E, 
the potential field varied more smoothly in the spanwise direction than did the poten- 
tial field obtained from the design solution. In light of the previous discussion, this 
result verifies that the inverse boundary condition was. in fact, forcing the adjacent 
potential fields away from each other. 

It should be noted that this problem is not due solely to the implementation 
of the direct-inverse technique since this oscillation has not been observed with the 
ZEBRA design code. Rather, it seems to be unique to the coupling of the method 
with the analysis code, FLO-30. Seemingly, tw r o pertinent differences between the 
two codes exit. Firstly, the ZEBRA code, w’hich uses a sheared Cartesian coordinate 
system aligned with the wing, applies the boundary conditions at the mean plane of 

2 , 7*7 
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Fig. 19 The potential field on the wing's surface for a diverging design solution 
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the wing. This first difference is important, since the actual thickness of the wing 
may have less of an impact on the flowfield computed by the ZEBRA scheme due to 
the fact that the point of application of the boundary condition is not changing with 
time. Secondly, its full potential, fully conservative numerical scheme uses a mid* 
segment type of finite difference approach rather than a finite- volume scheme with 
fourth derivative type compensation terms 26 that seem to be amplifying the errors in 
the design solution. 

Nevertheless, after exploring many alternatives to counter this oscillation prob- 
lem, four methods based on the previous observations have been devised to damp out 
the spanwise oscillation: 

A) Specify the inverse boundary condition at at least every other spanwise 
station and linearly interpolate the inverse displacements to the stations lying in 
between. This has been named the Type 11-2 method. 

B) Specify the inverse boundary condition at every station, but again only 
calculate inverse changes at every other station and linearly interpolate the inverse 
changes to the stations in between. This will be referred to as the Type II method. 

C) Immediately prior to every surface update, calculate all spanwise derivatives 
of the potential in the residual based upon a potential function smoothed in the 
spanwise direction. This smoothing is accomplished by first defining the operator <j^ 
as 

c d = \h - 1 + (i - I) fk + \h+ 1, 


0 < e < 1 


(4-2) 
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where e determines the amount of smoothing. Then using a in the spanwise differ 
entiation of <j> with the maximum amount of smoothing (i.e., 5 = 1 ) 


do 

KJij. 


= <T(8qO 


(4-3) 


the smoothed spanwise derivative of o becomes 


. ( oi,j.t-r: ~ Qi'.j.t d- Qj.j.k+i ~ 

W h j.h+ \ ~ 4.0 


(4-4) 


D) Smooth the slopes, in the spanwise direction in the design region in the 
same manner as with method C. It should be noted that, as stated earlier, smoothin G 
the integrated slopes, i.e the inverse corrections, did not suppress the oscillation but 
only slowed the rate of divergence. 

Three different cases were studied in order to test the effectiveness of each 
method at suppressing the oscillation and in reproducing the knovn target section. 
All three cases used Lockheed Wing- A at a Mach number of .8 and at an angle of 
attack of 2°. The first case utilized a NACA 0012 airfoil as the initial section and the 
original supercritical wing sections accompanying Wing-A as the target section. The 
design region stretched from 30-70% semispan of the wing and began 5% aft of the 
leading edge and extended to the trailing edge. Since a medium grid (80x12x16) was 
employed, there was a constant ( grid station at every 10% semispan. Results are 

2 . 8 ' 



shown in Fig. 20 for the four different approaches. 

Although all four approaches worked well for this case, by using the RMS of the 
errors between the target section and the section designed as a measure of accuracy, 
methods A and C produced the best results for this case in the interior as well as 
at the edges of the design region. For the same number of flowfield iterations, the 
technique D produced the most unsatisfactory results when compared to the target 
sections. 

The effect of each approach on the residual at the trailing edge after 10 surface 
updates can be seen in Fig. 21. The discontinuities in the residual for method A 
is due to the fact that the inverse boundary condition is applied only at the 30. 50 
and 70% semispan locations. All four approaches have a characteristic jump in the 
residual at the first spanwise design station at 30 % semispan. This jump is probabl} 
due to the previously discussed spanwise mismatch problem with the potentials at 
the direct-inverse interface, which manifests itself in the compensation terms. The 
Type II method had the largest jump at this interface, while the Type II-2 method 
had the smallest jump. Notice that the spanwise distributions of the residual for the 
two smoothing approaches are quite similar in the design region. 

Since only small differences existed between the methods for the previous test 
case, a more severe test was conducted by designing an entire wing using NACA 0006 
sections as the initial airfoils and NACA 0012 sections as the targets. These sections 
were chosen due to the fact that most of the problems in the past were amplified 
by beginning with a thin section and targeting a thicker section. Furthermore, a 
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Fig. 20 Comparison of airfoil sections designed using the four spanwise oscillation 
remedies on a medium gria from 30% to /0/o semispan. 
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flic four span wise oscillation remedies after 10 global il orations 





lull wing design would reveal whether the accuracy of each method depended on the 
spanwise location of the wing. 

When an attempt was made to compute these cases, it uas disco\ered 
when using the smoothtng two smoothing techniques, (methods C and D). it was 
necessary to use zero order extrapolation ol the displacements from the adjacent 
grid station to the root section. The root section tended to lag in the convergence 
process in comparison to the rest of the grid stations. Tins behavior is possibly due 
to a slowly converging flowfield at the the wing-fuselage juncture. Since all of the 
sections started out too thin, this lagging of the root section forced the adjacent 
grid station to quickly become too thick, which led to divergence at the root in both 
cases. Zero order extrapolation of the nondimensionalized displacements forced the 
root section to converge at a rate which was more in compliance with the rest of the 
grid stations at the expense of degrading the accuracy of the root section. Since the 
root section has been successfully designed independently, presumably, this problem 


might be circumvented by simply allowing the flowfield solution to converge further 
before each relofting, although such a procedure would probably be a less efficient 

approach. 

Also, no smoothing of the potentials or the residuals was used at the tip. Since 
both the residual and the potentials are quickly varying in the spanwise direction in 
the tip region, smoothing leads to large errors in the residuals and hence the section 
shapes. In fact, better results can be obtained for the smoothed potential approach 
by using a zero order extrapolation of the normalized displacements from the grid 
station inboard of the tip to the tip. Overall though, the inboard sections of the wing 
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slowly became thicker, while those outboard responded more quickly, initially causing 
these outboard sections to actually become too thick. 

The resulting sectional shapes for the four different methods are compared in 
Fig. 22. As can be seen in the figure, method C works well when designing in the 
interior of the wing, but did not give satisfactory results at the tip of the wing where 
smoothing the quickly varying potential led to large errors in the section shapes. Since 
the residuals also varied quickly at the tip, the slopes at the tip were not smoothed 
with method D. Since there w*ere not any slopes defined at the fuselage ghost point 
location, (z. ky, 2), the slopes w’ere not smoothed at the root either. This method 
produced the most accurate results w*hile still managing to suppress the oscillation 
problem. In contrast, the Type II and Type II-2 methods worked well on the entire 
wing surface, and nothing special needed to be done at the root or tip. 

The same case u*as executed on the fine grid (160x24x32) to study any effect 
of grid size on the accuracy* and effectiveness of the methods. This grid allowed 21 
design stations which w*ere located a distance of 5% semispans from each other. Vv hen 
using the Ty*pe II and Ty r pe IT2 methods, the lagging of the root section actually 
forced the section located at 10% semispan, two grid stations outboard, to become 
too thick, which led again to a divergent solution. Thus, for this fine grid case, it was 
necessary* to use zero order extrapolation of the the normalized displacements from 
the adjacent station to the root w*hen using all four remedies. Cases which do not 
require such large changes in thickness at the root have not required this procedure 
using the Type II and Type IT2 methods. In addition, because of the aforementioned 
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problems with the smoothing approaches' at the tip, no smoothing was used at the 
tip section. 

The results for this case are shown in Fig. 23. For this case, the smoothing 
approaches yielded satisfactory sections on the region of the wing spanning from 
about 30% to 85%. Elsewhere, the sectional shapes vary significantly from the target 
section. Thus, the smoothing approaches work well when designing in the interior of 
the wing, but they do not. give satisfactory results near the root and tip of the wing. 

An objective measurement of the accuracy of the sections in relation to the 
target can be obtained using a coefficient of determination, r, defined as 



where is the standard deviation of the ordinates of the target section defined as 

2 

[ 1 (y* ~~ ymean) 1 ^4 _ 0) 

V [ 71 — 1 J 


and 

2 

(y* ~ ydtsign) I £4 _ 7) 

• v '~ [ n — 2 J 

is the deviation of the design from the target for the same i values. This quantit} 
varies from 0 to 1, one being perfect. 

Moreover, to further clarify which method produced the least amount of oscil- 
lation, the average error variation in the spanwise direction for each method should 
be compared. The spanwise variation of the coefficient of determination and average 
percent error are shown in Figs. 24-29. The Type II and Type II-2 methods 
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Fig. 27 The cited, of the four spanwi.se oscillation remedies on the average percent 
error between the target and the design sections for the medium grid, full- 



Spanwlse Geometry Station 

Fig. 28 The effect of the four spnnwise oscillation remedies on the coellicient of 
determination for I. lie line grill, full wing design 




produced the least amount of oscillation, while smoothing the potentials produced 
the most amount of oscillation in the error. 

There is still some doubt by this investigator whether only secondary aberrations 
have been observed and not the true, fundamental cause of the oscillation. In light 
of this, another effect that should be investigated is that of the addition of mass into 
the flowfield by the inverse boundary condition. Some other investigators" 3,43 have 
included a source correction in the far field and in the near field 43 . In this research, 
this source correction was neglected since this addition of mass would be driven to 
zero at convergence. But, its effect on the unconverged solution is not clear. In order 
to see if this had a significant effect on the solution, a quick, numerical experiment was 
performed in which the distance to the outer boundary was doubled. (See Fig. 30) 
Presumably, if the addition of mass was adversely affecting the boundary condition 
in that region for a given distance, it would have less of an effect if the distance were 
increased since the additional mass flux arriving at the boundary would be less and the 
outer boundarv boundary condition would be better satisfied. When this computation 
was completed, however, the solution seemed to be completely unaffected, diverging 
at the same point in the iteration history. This was only a simple attempt at proving 
that the sources on the wing were not the fundamental motivation for the oscillation. 
A thorough analysis must consider the effect of this mass addition on the downstream 
boundary and the near field. The downstream boundary could be stretched further 
downstream, and appropriate source correction terms, using the flux ejected from the 
inverse regions of the wing as the source strength, could be added to the reduced 


potential in the entire flowfield. 
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Another possible cause could be the assumption of flow tangency used in the 
residual expression in the integration of the flow tangency equation. When this as- 
sumption is made, not only are the fluxes reflected about the current wing boundary, 
but so are the compensation terms. This procedure in effect doubles the amplitude of 
the and the 6 iv type compensation terms. Since the flow is not generally tangent 
to the current shape when designing a new airfoil section, reflecting the compensation 
terms may be initially incorrect. An alternative formulation may be needed. 

In retrospect, a few comments about the advantages of each method in different 
design situations are warranted. For instance, methods C and D give the designer the 
most flexibility; the desired pressure distributions can be imposed at every spanwise 
grid station, and the section shapes corresponding to each grid station can be calcu- 
lated relatively independently of the adjacent stations. On the other hand, because 
of the interpolation required in the first two methods, the section shapes at 'odd' 
stations are directly dependent upon the shapes at ‘even’ stations: so although the 
designer loses a little flexibility, he gains a smoother spanwise distribution of section 
thicknesses in the spanwise direction. From a designer's standpoint of course, method 
A is the most restrictive of the four, but it yields the smoothest designs in the span- 
wise direction, and converges the quickest. Therefore, method A (i.e., the Type II-2 
method) would most probably be the best to use with wings of moderate to high 
aspect ratios. But, Method B (i.e., the Ty r pe II method) would most probably be 
necessary for wings with aspect ratios of the same order as Lockheed Wing-B. 


lo } 
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CHAPTER V 

RESULTS AND DISCUSSION 

Since the versatility of the method in designing multiple, overlapping 
regions of the wing has already been well demonstrated by Gaily 44,13 , most of the 
test cases presented, herein, were chosen instead to exhibit some of the constraints 
and limitations of the current inverse design procedure. The cases were chosen to 
reveal the approximate limits imposed on the aspect ratio and sweep of the ving, and 
the significance of grid skewness, viscous interaction, grid refinement, and the initial 
airfoil on the final airfoil section design. Some questions about the compatibility 
of Mach number and pressure distribution will be answered by designing a wing 
at one Mach number using pressures obtained from a wing analysis at a different 
Mach number. Finally, preliminary results will be presented for a partial wing design 
beginning aft of the leading edge and terminating forward of the trailing edge. 

V.l Boundary Laver and Wake Effects 

One of the objectives of this study was to determine the significance of various 

viscous effect in the design of transonic wings 9 . The wing chosen for this study was 
a typical transport type wing, Lockheed Wing- A. This wing has an aspect ratio of 
8.0, a leading edge sweep of 27°, a taper ratio of .41, a twist of 2.28° at the root and 
-2.04° at the tip, and 1.5° of dihedral. 


3^3 
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An input pressure distribution was obtained by analyzing Lockheed Wing- A us- 
ing full viscous effects ; these included boundary layer displacement thickness, vake 
thickness, and wake curvature. The flight Mach number of .8. angle of attack of two 
degrees, and Reynolds number, Re, of 25 million used in the analysis were thought to 
represent flight conditions for a typical, average-sized transport; and the distribution 
was considered to be typical of that which would be available to and desired b\ a 
designer. All computations were performed on a fine (160x24x32) grid. The resulting 
pressure distributions obtained from the analysis were used in two separate design 
cases, each composed of five and three subcases, respectively. The first series of cases 
was a full wing design using the target section as the initial section. By using the 
target section, any effect of the initial section on the final outcome would presum- 
ably be eliminated. The type II design method was used ano the inverse boundar} 
condition was enforced from 5% aft of the leading edge to the trailing edge. Further- 
more, relofting was not initially done at all. The results for the partiall} converged 
cases were plotted and then further converged allowing relofting to take place. In 
this way, the effect of relofting on the final design could be determined. The itera- 
tion history of each case was kept the same, even though by doing this the absolute 
level of convergence could ver)' well be different since changes of various magnitudes 
were associated with each case. The large amount of computational time required for 
these cases dictated this type procedure and for comparison purposes this approach 
is acceptable. Fortunatel} r , it turned out that the sectional shapes in every case were 
varying quite slowly by the end of the design run, indicating that the sections were 
near convergence. 
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In each case, viscous options were 'turned off’ one at a time to assess their 
effect. In the first case, the wing was designed with all viscous effects. In the second, 
the lag entrainment was turned off. The third case did not use wake curvature, while 
the fourth neglected both wake curvature and wake thickness. Finally, in the fifth 
case the wing was designed inviscidly. The resulting umelof.ed designs for each case 
are compared in Fig. 31. As expected, the inviscidly designed sections are slightly 
thicker at the root where the l.ormalited boundary layer displacements are thinnest 
(see Fig. 32 ) and become increasingly thicker towards the tip in accordanc 

the thickening boundary layer. 

Neglecting lag entrainment, wake curvature and thickness had very little effect 
on the designed sectional shapes overall. But. if the trailing edge region is examined 
closely for cases with the wake effects neglected, the trailing edges sometimes cross. 

Upon converging these shapes further and enforcing a trailing edge ordinate 
requirement with relofting, significantly different results were obtained. As shown 
in Fig. 33, the inviscidly designed shapes are now thinner on the upper surface and 
slightly thicker on the lower surface, especially in the cove region where v iscous e..ects 
are large. Also, because of the relofting involved, the leading edge radius has become 
smaller. The rest of the cases produced sections which did not deviate much from 
target, except near the tip. However, neglecting both wake effects produced sections 
that were actually thicker than the target. This change was due to the relofting that 
was necessary to uncross the trailing edges, which produced larger leading edge radii 

and hence thicker sections. 

3 
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Fig. 31 Comparison of the unrelofted sectional shapes designed using different vis- 
cous interaction assumptions with a pressure distribution obtaine^d from 
a fully viscous analysis of Lockheed Wing- A at a Re = 24 x 10 . M — 
.8, and an a = 2° 
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Table 1. 


Results from the analysis of the wings designed with different viscous in- 
teraction assumptions at a M = -8 and a Rc = 24 x 10 6 


Case 

Wing C'i 

Wing 

— Fuselage C'i 

Cd 

Target 

.4745 

.5347 

.0197 

Full Viscous 

.4636 ! .5226 i .0195 

! \ 

No Lag Entrainment 

.4719 

.5316 

.0197 

No Wake Curvature 

.4636 

.5226 

.0195 

No Wake Effects 

.4605 

.5194 

.0193 

Inviscid 

.4060 

.4598 

.0169 


The resulting wing for each case was analyzed using full viscous effects and the 
same iteration history. Table 1 gives a comparison of the lift and cirag coefficients 
resulting from the analyses of these designed wings. 

As can be seen from the pressure distributions shown in Fig. 34 and Table 1, 
the mviscidly designed wing produced 15% less lift than did the target wing. The lift 
usually obtained in the cove region was diminished, in this case, by the decambering of 
the aft portion of the wing. The thinning of the top in conjunction with the thickening 
of the bottom of the inviscidly designed airfoils also caused a decambering of each 
section, which explains the large decrement in lift produced. As shown in Fig. 35, the 
reason the top was thinner is because the boundary layer displacement thicknesses 
which are ‘built’ into the imposed pressure distribution were not subtracted from 
the inverse displacements in the inviscid design. In order to meet the trailing edge 
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ordinate requirement, the resulting section had to be relofted more to compensate, 
thus leading to a thinner section on top. 

The wing lift coefficients obtained from the analyses indicate that by not using 
lag entrainment, a design correlating closely with the target can be better accom- 
plished for the given sequence and number of flowfield iterations. It is suspected 
displacements and hence the inverse displacements may take longer to converge to 
the correct value as compared to excluding lag entrainment. By ignoring wake curva- 
ture and using all the other available viscous options, wings with identically slightly 
lower lift coefficients as compared to the targets were produced. Furthermore, wake 
thickness influenced the design in a slightly more profound way than did wake curva- 
ture bv producing a wing with 3% less lift than the target. 

As an after thought, the original wing was analyzed with each viscous option to 
assess its effect. The analysis results of the designed wings, shown in Fig. 36, reveal 
that wake curvature effects were practically negligible. This result may be due to the 
relatively high freestream Reynolds number of 25 million used in the comparisons. 
Since this Re would lead to low values of 6 ' and 0, the curvature effects would also be 
expected to be low; Streett's case 20 used a much lower Reynolds number of 6 million. 
On the other hand, neglecting wake thickness and lag entrainment effects both had a 
decremental effect on the wing’s lift, which was probably due to the forward shifting 
of the shock location. 

The second set of design cases involved a partial wing design which extended 
from 30-70% semispan and began 10% aft of the leading edge ot the airfoil, but the 

3 1 L 1 
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inverse boundary condition was only enforced at the 30, 50 and 70% semispan station 
and the displacements were linearly interpolated to the stations in between. The 
initial airfoil section at 50% semispan was formed by thinning the supercritical target 
section by 6% and removing the cove region. The initial sections at the edges of the 
design region were the same as the target sections, while the remaining sections were 
obtained through linear interpolation. The results for these cases are presented in 
Fig. 37. For the Reynold’s number chosen, neglecting wake effects seems to ha\e 
had a small effect on the resulting design. The sections are a little thicker than the 
sections designed with full viscous effects. As noted earlier, the wake effects had 
relatively little effect on the pressure distributions obtained from the analysis of the 
target wing: but, when the boundary layer displacement thicknesses' obtained were 
investigated, it was discovered that neglecting wake effects in the analysis produced 
boundary layer displacement thicknesses that were on the average 3.5% thicker at the 
trailing edge than those obtained from a full viscous analysis. Since the boundary 
layer displacement thicknesses are subtracted from the initial inverse changes to }ield 
the hard airfoil, these larger displacement thicknesses would produce a section that 
was initially thinner than the target; but, after relofting the airfoil section, it would 
actually be thicker than the target. 

The wing sections designed inviscidly are profoundly different at 30 and /0% 
semispan, but only slightly different at 50% semispan. The thinning of the top surface 
in complement with the thickening of the lower surface significantly decambered these 
sections. The large differences at the inboard and outboard design stations are due to 
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Fig. 37 Comparison of the relofted sectional shapes designed using different viscous 
interaction assumptions with a pressure distribution obtained from a fully 
viscous analysis of Lockheed Wing- A at a Re = 24 x 10 6 , M = .8, and an 
a = 2°. The design region extended from 30%-70% semispan. 
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ihe influence of the inviscid pressures outside the design region; and, the remarkable 
agreement in the middle of the design region, except in the cove region where the 
boundary layer is thick, is due to the influence of the viscous boundary condition 
at the edges of the design region. This observation can be verified by reviewing the 
previous case and noticing that the airfoils sections varied smoothly in the spanw.se 
direction M all spanwise stations. 

After the wings were designed, all three were then analysed with full viscous 
effects to assess the significance of the changes made to the wing on the pressure distri- 
butions and to see how well these pressures matched the target pressures. Knowing 
that the wing designed with full viscous effects is correct, it is quite obvious from 
Fig. 38 and Table 2 that the wing designed inviscidly is quite unsatisfactory. The 
shock is not far enough aft and the lift produced is sometimes 20% smaller than that 

desired. 

Based on the results of this study, it can be concluded that for the Reynold s 
number and Mach number chosen, wake curvature and wake thickness and lag en- 
trainment have a very small effect on the designed airfoil sections. However, the 
boundary layer displacement effect has a profound effect on the section shapes and 
hence must be included in the design process to yield a wing which will produce the 
desired lift in a viscous environment. 


V.2 Span vise Grid Skevness 

In the course of the present research, it vas discovered that the skevness of the 
grid lines leaving the tip of the ving (Fig. 39) can have a dramatic efiec 
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Fig. 38 Comparison of pressure distributions obtained from a viscous analysis of 
the Lockheed Wing- A using three different viscous interaction assumptions 
with a M = .8, ct = 2°, Re = 24 x 10 6 
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the design of the sections near the wing tip. As can be seen in Fig. 40, if the grid « as 
significantly skewed and the input pressures were calculated on an nonskewed grid, 
i, was impossible .0 obtain the correct airfoil shapes in the tip region. This difficult 
is due to the large differences in pressures between the skewed and nonskewed grid. 
These pressure profile differences are shown in Fig. 41. As shown in the fi 0 ure. 
the grid skewness has caused the shock location to move further aft. Although the 
skewness of the grid was quite extreme in this case, these results affirm the need for 
smoothly varying grids in wing design, at least in the spanw.se direction. It should 
be noted though, that if the input pressures were obtained on a skewed grid and used 
in the design process with a skewed grid then the tip sections were well resolved. In 
summary then, if the pressures calculated on an nonskewed grid are correct or closer 
to real pressures encountered in flight, then it would be wise to ensure that the gnd 

is smoothly varying. 


V.3 Wing Planform Effects 

Three cases were attempted to roughly delimit the applicable range of aspect 
ratios and leading edge sweep angles for which good results could be obtained with 
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Fig. 40 Sections designed with a skewed grid using pressures obtained from an 
analysis on a nonskewed grid 
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Fig. 41 Comparison of a pressure distribution obtained from an analysis of a skewed 
grid with one obtained from an analysis of a nonskewed gn 
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Fig. 42 Grid generated about Wing-C with an incompatible root section and fuse- 
lage cross section 

the present design method. These included Lockheed Wings A. B and C. These wings 
have aspect ratios of 8. 3.8, and 2.6. leading edge sweep angles of 27, 35, and 45 degrees 
and taper ratios of .4, .4, and .3 respectively. The target pressure distributions were 
obtained by a direct analysis of the target wings in an inviscid environment. The 
initial section for Wing- A was a NACA 0012, while a NACA 0006 was used for Wing- 
B. The original section was used with W : ing-C due to the difficulty of the case. Also 
for Wing-C, as opposed to the circular cross-section, an elliptical cross section of the 
fuselage was used to provide a flatter surface for the grid generation package. The 
circular cross-section combined with the large relative thickness of the root section 
compared with the width of the fuselage played havoc on the grid at the root, as can 

3 H 


be seen in Fig. 42 
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1„ order to better understand the flow about each wing, the corresponding 
velocity vectors on the surface of each wing were plotted, as shown in Figs. 43-45. 

As should be expected, the spanwise component of the flow increases as the aspect 
ratio decreases and sweep increases. It is also interest.ng that there seems to be 
an inboard component of the flow for all three cases on the upper surfaces aft of the 
leading edge. Tins .aboard flow may be attributed to the effect of the fuselage and the 
wing tip vortex. These effects can be seen most readily by viewing a cross section of 
the flow just aft of the wing tip shown in Fig. 46. The vortex near the tip of the wing 
is quite evident, and flow langency at the fuselage also contributes to the spanwise 
component of the flow. The momentum of the air over the tip must dominate the 
flow, since, as seen in Figs. 47-49, the spanwise pressure gradients appear to encourage 
the air to move outboard. However, in order to determine whether 

the flow actually traveled in the inboard direction, it would be necessary to plot the 


actual streamlines of the flow over the surface of the wing. 

The design region for Wing-A and Wing-B extended from 10-100% semispan 
and began 5% and 2.5% aft of the leading edge, respectively. Computations were 
performed on a fine grid. Results for Wing-A are shown in Fig. 50, while results for 
Wing-B are shown in Fig. 51. As can be seen the designed and target sections for 
both wings are in excellent agreement in the interior of the design region and closely 

match at the edges of the design region. 

In the case of Wing-C, the section shapes should not have changed with the 

application of the inverse boundary condition. But, because of the large amount 
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Fig. 46 Cross-section of the flow in the x-y plane for Lockheed Wing-C 
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Fig. 48 Pressure contour plot for Lockheed Wing B M 
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Fig. 49 Pressure contour plot for LockWd Wing-C M 
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Fig. 50 Comparison of the designed sections with the targets and the initial sections 
for a fine grid case using Lockheed Wing-A 
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Fig. 50 Continued 
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of spanwise flow and the associated spanwise gradients for Wing-C, the spanwise 
oscillation effect could not be overcome with an)' of the present remedies. Further 
information about this case was obtained by using the Tvpe II method and not 
relofting the section shapes. The results for such a converging fine grid case are 
shown in Fig. 52. The first design station at 18% semispan is too thick on 

the upper surface as compared to the target. This discrepancy is again due to the 
o\ er prediction of the residual at the first station due to the initial mismatch in 
the potentials in the spanwise direction, and. hence, to large spanwise gradients of 
the potential. The errors diminish as the tip is approached, but are always relatively 
large in the trailing edge region due to the difficulty in accurately imposing the inverse 
boundary condition near the trailing edge for this case. If an attempt were made to 
converge this case further by continuously reloft-ing the shapes to meet the trailing 
edge ordinate, the same spanwise oscillation problem would again occur. However, 

non- relofted results such as in Fig. 52 would be very useful for preliminary design 
studies. 

V.4 Initial Profile Effects 

One of the disadvantages of the direct-inverse method is that a priori knowledge 
about the correct shape of the leading edge must be known to achieve suitable airfoil 
shapes and desired trailing edge thickness. Relofting does alleviate this disadvantage 
to a large degree; but it will not, in general, produce a leading edge that will yield the 
desired pressure distribution at the leading edge if the inverse boundary condition is 
by necessity applied too far aft. It was thought that because FLO-30’s grid package 
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Fig. 52 Comparison of the designed section with the target for an unrelofted fine 
grid case using Lockheed Wing-C 
















clusters grid lmes close to the leading edge of the airfoil, that the design could be 
started quite clo>e to the leading edge, thus relieving the designer of the difficult}* 
of choosing a correct nose shape. Two test cases were conducted to investigate the 
dependence of tie final design on the initial airfoil section. Both used Lockheed 
Wing-A at the some conditions mentioned earlier for the viscous study. For the first 
case, the initial airfoils were the same as those in the viscous study. These airfoils all 
had leading edges which were in the same family as the target section. The design 
began 10% aft of the leading edge. In the second case, NACA 0012 sections were 
used at all the .design stations; here, the leading edge of these sections were not in 
the same family as the target airfoil sections. For this case, the pressure boundary 
condition began 4% aft of the leading edge. Referring to Fig. 53, it can be seen that 
although slightly better results were obtained near the leading edge for the first case, 
that the airfoils designed were fairly insensitive to the initial section. 

V.4.1 Direct-inverse interface proximity to leading edge 

Since experience with the method has shown that the closer the inverse bound- 
ary condition is applied to the leading edge, the longer it takes for the solution to 
converge, it was of interest to determine how the location of the direct-inverse in- 
terface affected the final design and the resulting pressure distributions. This study 
was accomplished with the aid of the previously discussed Wing-B case, whose design 
region began at 2.5% chord, and an inviscid design of Wing-B also with NACA 0006 
sections as the initial geometry. With the second case, the design was started at 5% 
chord from the leading edge; and, the input pressures were obtained from an inviscid 
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Fig. 53 Comparison of sections designed using two different initial sections 
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analysis of Wing-B. Since the design pressure distributions were consistent in both of 
these cases, the fact that one was a viscous design and the other an inviscid design is 
not important here. 

Some representative samples of the resulting section shapes for the second case 
are shown in Fig. 54. The resulting wings were analyzed under the same conditions 
that the original input pressure distributions were obtained. Representative samples 
of the resulting pressure distributions are compared to their respective target distri- 
butions in Figs. 55,56. As can be seen, the wing whose design began 2.59c aft 
captured the suction peek at the leading edge, while the other case, which began at 
b% aft of the leading edge, did not. 

When designing near (less than 59c) the leading edge, the solution sometimes 
began to slightly diverge or ceased converging. Lsually the design could be converged 
to the point where there was only a maximum change in the surface of chord. 

This was more a problem on the fine grid than on the medium. If it was necessary to 
converge it further, the beginning of the design region was moved aft. This observation 
is important because if it is necessary to begin the design close to the leading edge to 
properly determine the shape of the nose, a successful design may be accomplished by 
beginning the design as close to the leading edge as desired or is possible, then moving 
the beginning of the design region aft as the solution approaches the last stages of 
convergence. This method not only frees the designer from the task of choosing the 
correct leading edge shape, but it should also accelerate the convergence of the design 
considerably. 
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Fig. 54 Comparison of the designed sections with the targets and the initial sections 
for a fine grid case using Lockheed Wing-B and a design region beginning 
at 5.0% aft of the leading edge. 
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Fig. 56 Comparison of the pressure distributions obtained from an analysis of the 
Wing-B design, which had a design region that began 5.0% aft of the 
leading edge, with the target pressure distributions 
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Because of the leading edge clustering of grid points in TAW5D, successful 
designs have been accomplished on the medium grid with the chord wise direct-inverse 
junction beginning just aft of the stagnation point on the lower surface. If the pressure 
boundary condition is applied upstream of the stagnation point, major difficulties 
arise when an attempt is made to integrate past this point of singularity, since the 
slope, jt, is indeterminate there. 

For the case shown in Fig. 57 . the design was begun 1% aft of the leading 
edge, but in retrospect, it could have begun close to .3% aft of the leading edge since 
the converged stagnation point was located about .29c aft. Notice how precisely the 
designed surfaces can be computed when compared to the targets outboard of the 
first design station. This case effectively demonstrates that since the design region 
can be extended extremely close to the leading edge with TA^ 5D. the fact that the 
pressure boundary condition can only be applied aft of the leading edge is a \er} 
small shortcoming of this direct-inverse method. 

V.5 Pressure Distribution Compatibility 

Since a designer might not readily have available an input pressure distribution 
compatible with the design freestream Mach number, the effect of designing a wing 
at one Mach number using a pressure distribution obtained from an analysis of the 
wing at a different Mach number was investigated. The Wing-A planform was used 
throughout this portion of the study. NACA 0012 sections were used as the targets 
and NACA 0006 sections were used as the initial sections in the design. The entire 
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Fir 57 Comparison of the designed sectrons with the targets for a case whose 
° design regton began 1% aft of the leading edge 
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wing was designed on from root to tip, and the design region started 10% aft of the 
leading edge of the wing. 

Two separate tests were performed. The first involved a fine design at a nearly 
incompressible Mach number of .2 using a pressure distribution obtained from an 
analysis of the target at a Mach number of .1. As can be seen from Fig. 5S. 
thinner section shapes were obtained at the higher Mach number. This thinning is in 
agreement with the 2-D Prandtl-Glauert similarity rule 45 

:= = £5 , 5 -D 

T2 v /l - Ml 

which states that the C p will be invariant with Mach number if the thickness. r. is 
reduced as the Mach number is increased for linearized flow. For this case. Eq. (5- 
1) would predict that a 1.54% decrease in thickness would be necessary to have the 
same pressure distribution at the higher Mach number. The design code for this 3-D 
case produced a section which was on the average 1.6% thinner than the NACA 0012 
section. 

The second case involved a medium grid design at a Mach number of .85 using 
a pressure distribution obtained at a Mach number of .80. Referring to Fig. 59, the 
section shapes produced are again thinner than the initial section. The top surface, 
though, required a sudden thinning of the surface at the shock location. Surprisingly, 
upon analyzing this wing, the pressure distributions shown in Fig. 60 match quite 
well with the target everywhere except in the tip region of the wing. So, given the 
constraints of the problems, it appears that the only way the boundary conditions 
could be met was to have these dips in the airfoil surface. Since these dips might 
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Fig. 58 Comparison of the section designed at a M - 
distributions obtained from an analysis at a M 
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Fig. 58 Continued 
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Fig. 59 Comparison of the sections designed at a M - .85, using input pressure 
distributions obtained from an analysis of Lockheed "Wing- A with NACA 
0012 airfoils at a M = .8, with the original NACA 0012 sections and initial 
NACA 0006 sections 
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Fig. 60 Comparison of the pressure distributions obtained from an analysis at a 
M = .85 of the Lockheed Wing-A which was designed at a M = .85 using 
input pressure distributions obtained from an analysis at a M = .8, with 
the input pressure distributions 
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lead to boundary layer difficulties, it would probably behoove the designer to vary 
the Mach number or alter the pressure distribution to eliminate the necessity of these 
dips. 

V.6 Grid Refinement Effects 

Since the computational time required for a design on the medium grid is about 
an eighth of that required on a fine grid, it may be tempting to try to design on the 
medium grid using fine grid or real pressures. In order to assess the practicalitj of this 
approach, a transonic design on a medium grid using fine grid pressures was carried 
out. The case was performed at a Mach number of .8 and an angle of attack of two 
degrees. The original supercritical sections for Wing- A were used as the initial, as 
well as. the target sections. The results are shown in Fig. 61. The only place vhere 
the designs came close to the target was near the middle of the wing. A slight va^ 
appears in the upper surfaces of the designed sections near the shock location. This 
pertubation is due to the smearing of the shock on the medium grid. The section 
designed at the wing tip deviated considerably from the target. The fact that at the 
wing tip the fine grid C\ is lower than the medium grid C\ most probably led to the 
decambering of the sections at the wing tip. 

No attempt was made to match the Ci's of the fine grid and medium grid 
analyses by varying the Mach number or angle of attack, but a comparison of the 
medium grid pressures at various Mach numbers and angles of attack with the target 
fine grid pressures for the supercritical wing shown in Fig. 62 reveal that it would 
probably be necessary to alter the twist of the wing to closely match the Ci s at all 
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of the design stations. It also shows that increasing the angle of attack to 2.1° would 
have produced closer matching Cfs and hence perhaps better designs. In retrospect, 
though, given that the fine grid pressures are correct or more realistic, it would be 
necessary, unless appropriate corrections can be found, to use the fine grid to proper]} 
design the correct airfoil sections. 

Y.7 Fixed Trailing Edge Design 

This case was investigated to verify that a fixed trailing edge design could be 
accomplished with the present version of the code. The case chosen utilized Lockheed 
Wing- A at a Mach number of .8 and an angle of attack of 2°. A NACA 0012 section 
was used as the initial geometry from 30% to 70% semispan, while the remaining part 
of the wing used the original supercritical sections. The inverse boundary condition 
was enforced from 5% to 80 % chord. The airfoil aft of 80% chord was fixed so that it 
maintained the NACA 0012 trailing edge shape. The input pressures were obtained 
through a medium grid inviscid analysis of the wing with the original supercritical 
sections used throughout. Furthermore, to provide for a smooth transition at the aft 
direct-inverse junction, the displacements were smoothed in the chordwise direction. 
The type II-2 design method was used in this case. 

The resulting section shapes are shown in Fig. 63. The target airfoil section 
■would actually be the first 80% of the supercritical section and the last 20% of the 
NACA 0012 section. Surprisingly, even with the aft portion of the wing fixed, the 
designed sections came quite close to matching the original Wing- A profiles at the 
30% and 50% semispam locations. At the 70% semispan location, the designed section 
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Fig. 63 Comparison of the sectional shapes designed with a fixed trailing edge 
region with a NACA 0012 and the original Lockheed Wing- A section 
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as compared to the original Wing- A section is much thicker on top and thinner on 
the bottom leading to a more cambered profile. This shape is probably due to the 
interaction of the geometric constraints and the required design pressures. The shock 
strength of the input C p distribution does become quite large at this location and 
it appears that the section may have become more cambered to account for this 
increase. Or. the increased camber may have been needed to provide the necessan 
lift required bv the inverse boundary condition. The pressure distributions obtained 
from an inviscid analysis of the resulting shapes are compared with those produced 
bv the original Wing- A sections and the NACA 0012 sections in Fig. 64 The figure 
reveals that the design pressure distributions are a combination of the Y\ mg- A ana 
NACA 0012 pressure distributions. It is also interesting that it seems a secondary 
shock near the aft limit of the design region was necessary to meet the constraints of 
this problem. This very impractical case, of course, was only meant to demonstrate 
that it is feasible to fix the aft region of the wing. If a more realistic trailing edge 
were used, better results would surely follow. 
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CHAPTER VI 

CONCLUSIONS AND RECOMMENDATIONS 

Progress in the direct-inverse wing des.gn method in curvilinear coordinates has 
been made, which included the remedying of a spanwise oscillat.on problem and the 
assessment of grid skewness, viscous interaction, grid refinement and the initial airfotl 
section on the final design. Some of the important conclusions were: 

(1) In response to the spanwise oscillation problem, designing at every other span- 
wise station produced the smoothest results for the cases presented. 

(2) A smoothly varying grid is especially needed for the accurate design at the wing 

tip. 

(3) The final designed airfoil section is independent of the initial section 
chordwise direct-inverse junction is moved close to the leading edge. 

(4) Boundary layer displacement thicknesses must be included m the successful 

design of a wing in a viscous environment. 

(5) Presently the design of only high and medium aspect ratio wings is possible 

with this code. 

(6) A partial wing design beginning aft of the leading edge and terminating p 
to the trailing edge is possible with the present method 

(7) Designs must be performed on a fine grid. 

It is recommended that more work be done to fully understand the fundamental 
motivation behind the spanwise decoupling problem in order to eliminate all spanvuse 
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oscillations in sectional thickness from the solution. This work should also include the 
development of a better way to handle the formulation of the residual at the span wise 
direct -inverse junction to eliminate the initial spanwise jump in the residual located 
there. Furthermore, the design scheme at the wing root and tip should be refined to 
provide more accurate airfoil sections in those regions. 

In addition, the necessary logic should be added to begin the integration of the 
flow tangency boundary condition on either side of the section s stagnation point at 
the present iteration level. This addition should allow the entire airfoil section to be 
designed with the pressure boundary condition specified everywhere on the ving s 
surface except at the stagnation point. 

Prelimenary results have indicated that by allowing the trailing edge ordinate 
to float an untwisted wing can be twisted. If this is a well-posed problem, methods 
should be devised to accurately calculate the twist given the inverse displacements at 
the present time level and to include this in the iterative process such that the twist 
angle converges without undue oscillation. It would also be interesting to investigate 
the possibility of also allowing the leading edge ordinate to vary m a constrained 

fashion so that the local dihedral angle could change. 

And finally, since the potential solution and, hence, the design, converge rather 
slowly due to the SLOR numerical scheme, the design scheme should be incorporated 
into the multi-grid version of FLO-30 to hasten convergence. 
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APPENDIX A 

DERIVATION OF THE FULL POTENTIAL 
EQUATION IN CURVILINEAR COORDINATES 


The full potential equation transformed from cartesian to curvilinear coordi- 
nates is derived here as a courtesy to the reader. 

The full potential or the continuity equation written in cartesian coordinates is 


(p u )z + (P v ) v + {pw)- - 0 


(. 4 - 1 ) 


where 



It is desired to transform this equation to a curvilinear coordinate system of 


77 , and ( where 


£ = £(*»?,-) V = v(z,y-.z) ( = ((x,y,z) 


(.4 - 3 ) 


By using the standard chain rule, the following operators can be defined 

d , d Od 

r = £=7 + vz - - t ( z 7 

1 ( v c 

d e d d d 

y -tv{- r7 lv-+<V( (- 4 “4) 

d _ . d , d a 

7 ~ ~ Vz~ + (ij 

2 i v C 

Using these operators in Eq. (A-l) yields 

(= (p u )( + Vz (pu) v + Cz ( pu) ( 

+ (y (p v )( 4- Vy {p v ) v -1- Cy (P v )( (.4-5) 

-f i s ( pw ) i -f Tj z (pu<) v -r (z (pw) ( = 0 
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Defining the Jacobian, J , as 


J = 




(t (y 
Vz Vy Vz 
Cz Cv Cx 


(*-6) 


d[x, y,z) 

Then after Holst, multiplying the Eq. (A-5) by J~\ and rearranging to conservative 
form plus remainder gives 

[((pu)^J- 1 ) + ((/n'KV 1 ) + ((^)^- 7_1 )]f 
4- [((pu)^J- 1 ) 4- ((^)V' 1 ) + ({p™)V:J-')} v 
4- [((pix)CxJ -3 ) + ((H Cy^" ] ) + ((^)C--^ _1 )] ( 

- (py) + (vzJ~ : ) v + {(=J *)<] 

- (p v ) (&-')( + (Vy^) v + (tyJ 3 ) ( ] 


- ( P w) + b*J’\+ M _1 ) c 


= 0 


Now using the fact that 


SJ- 1 _ 7 oJ_ 

ds ~ 6s 


(-4-7) 


(.4-8) 


the last three terms in brackets can be shown to be zero. For example, equating the 


first of these terms to zero 


(pv) (C= J *){ + (jlzJ l ) v + i(z J ^c ]- 0 


and expanding the derivatives and collecting like terms gives 

J 1 (£s){ + {Vz) v + (Cz)(] 


- J ~ 2 [tzJ( + VzJ v + <^c] = 0 


Rewriting Eq. (A-4) in matrix notation 


f 


fu 

Vs 

C. \ 

(1 

\ 

d 


= Cv 

Vy 

Cv 

V 

V 


u 

1 

\iz 

Vz 

Cr J 

(? 

J 


5 81 


(.4-9) 


(A -10) 


(A -11) 



After solving for Jji 5^, this becomes 


(n \ 

d_ 

OT ) 

U ) 

where : 


An -An An 

—Ai 1 >12; —A23 

A 3 \ — A 32 A 33 



(.4-12) 


-4n =Vy(: — CyVz An = (y(t — (z(y An = (yVz — T)y(z 

An =Vz(z - V:Cs -422 = (zCt - ( : ( z A 23 = (rf. - ( z 7) z (.4 — 13) 

-43) =Vz(y - Vy(>: A 32 = £zCy ~ ty(z An = ( z T]y - 7),( y 

These operators can be used to expand the derivatives of jj_, and (. so that 


Hz) l — [Ani- Z — Azi^-y -r Ani xz ) J 1 


(Vz) v - [AnVzz ~ A 22 Tj~y -r A 32 t]zz}J 1 • (.4 - 14) 

(C=)( = [.4)3C=; - -423C=v 4- AzzCzzjJ- 1 
Substituting these into Eq. (A-10) and collecting terms jdelds 

J '[Anizz — -4:)^ =y -f A 3 \(zz 
-4l2^z: — -422^=y T -432s=r 


-r-4)3^== - -423^ry T -433^-] 

-J~ 2 {An JzU - A 2 \Jyiz + A 31 JzU 
—AnJ=Vz — A 22 J v rj. -r A 32 J z t). 

—AnJxCs — A 23 JyCz "T -4 3 3JrC=] = 0 

with J = (.A n - TjzAn + ( X A U 

Expanding the second term in brackets in the previous equation to 

— J 3 [J Z J T J y ((~T] z (z. — (xVzCz 

+t*V*(z - (zVzCz -f (zVzCz - UVzCz) 

—JzitzVzCy — CzVyCz — ZzVzCy 


(A - 15) 


(.4 - 16) 


HyVzCz -r (zVyCz ~ tyVzCz)] 


1 



and cancelling like terms, this reduces simply to 


Second term — J 2 Jx 


where : J x =(* (An)j — ’J* (^12)* + Cx (^13)1 
+ txzM\ ~ VzzMl + Cxt-^13 

Partially expanding J z to 

J z = ^rr-4n — ’JxxAn ~ C=r-^13 

-viz (VzyCz + VyCzz ~ Vztzy ~ (yVzz) 

-Iz {izyCz -T iyCz: ~ tzzty ~ CrCxy) 

-r( z [tzyVz + (y r lzz — Vzytz ~ Vytzz) 

Upon collection of like terms, this becomes identical to the first term ii 
Eq. (A-15), thus satisfying the equality. This can be shown to be true 

remainder type terms in Eq. (A- 7). 

I'Cow, reducing the conservative part of Eq. (A 1 ) to 
[J -1 ((/m) U + ipv) ty + {pro) (z)] f 

+ [j- 1 ((/ni)i? s + (H J ?» + (/ w )^)], 

+ [J - 1 ((/Hi) C= + (pv) Cv + (p™) &)] c = 0 
and defining the contravariant velocities, U,V,W, as 


U \ (U tv tz 

V = [ Vz Vy Vz 

W J \Cx Cv C* 


u 

V 

. w 


with 


z ( z v z ( 

h = J~ l = V( Vn VC 

. z i z n z c 

Eq. (A-20) reduces to the desired conservative form of 


(phU), A (phV), + (phW) ( = 0 


(A -17) 
(A -18) 

(A -19) 

t brackets in 
of the other 

(A - 20) 

(A -21) 
(A - 22) 

(A - 23) 



APPENDIX B 


DERIVATION OF THE C p EQUATION 


Although this derivation probably appears in most good books on aerody- 
namics, it is included here as a courtesy to the reader. 


C T is defined as 


C p — j 


P- Poo 


. , (5-1) 

2 Poo? oo 

Using the definition of the speed of sound and isentropic relations, this can be rewrit- 


ten as 


Cp ~ ~ f M- 


F DO 


- 1 


(5-2) 


It is desired to obtain a relation for the pressure coefficient, C P , in terms of soley 
the freestream Mach number and the local q^. This can be easily accomplished by 
beginning with Eq. (2-14), 

p = (aM 00 )^ (B — 3) 


and using the isentropic relation in Eq. (2-10), pressure can be written ?.s 

- = (5-4) 

Pod 

Upon substituting this into Eq. (B-2), equation, Cp becomes 

Cp = 7 ” l ) ( B ~ 5 ) 

And finally, making use of Eqs. (2-8) and (2-9), the previous equation can be reduced 

to the desired relation : 


y 



(B — 6) 


.htK 5 J = (, u 2 -f v- + W-) . 
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Appendix IV 

User's Manual for TAW5D 


Originally, it was planned to include the user's manual 
•final report. However, due difficulties in the preparation of 
it will not be available until mid-November. Thus, it will be 
to NASA Langley under separate cover. 


in this 
f i gures , 
prov i ded 
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